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Preface 


The purpose of this SpringerBrief is to review heat transfer in layers of convective fluid. 
Six different configurations are considered—three that are versions of Rayleigh-Benard 
(RB) convection, which is driven by differential heating at the boundaries, and three 
that are driven by uniform internal heating. The essential features of all six models 
are derived mathematically. The experimental literature is reviewed in depth for the 
models of internally heated (IH) convection, which are much less studied than their 
RB counterparts. Experiments on RB convection are treated in less depth, as they have 
been thoroughly reviewed elsewhere. 

Along with placing the various convective models within a conceptual framework 
that brings out their similarities, we give some minor results not published elsewhere. 
For instance, a few of the linear instability and energy stability thresholds given in 
tables 2.2 and 2.3 either have not been reported before or have been reported with less 
precision. One of the bounds proven in §2.3| is also new, as are the visualizations of 
simulations included as figure [L^ 

Chapter[T]provides background and then defines the six configurations under study, 
the governing equations of our models, and their basic features. Chapter presents 
results that can be derived mathematically from the governing equations: linear and 
nonlinear stability thresholds of static states, along with proven bounds on mean tem¬ 
peratures and heat fluxes. For the IH cases only, chaptergives a quantitative survey 
of heat transport in both laboratory experiments and numerical simulations, followed 
by suggestions for future work. 

The author is grateful to Charles R. Doering, Erwin R van der Poel, Jared P. White- 
head, and Francis A. Kulacki for their many helpful comments on the manuscript. 
Thanks also to Martin Worner for providing his original data and Francis A. Kulacki 
for providing not only his original data but also many hard-to-find references. More 
general thanks are due to Edward A. Spiegel, who taught the author much of what he 
knows about convection and many other topics. 


David Goluskin 
University of Michigan 
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This review is concerned entirely with convection—fluid motion driven by differen¬ 
tial body forces. We focus on several simple configurations that lend themselves to the¬ 
oretical and experimental study. Convection arises in many contexts and figures promi¬ 
nently in astrophysics, geophysics, and certain engineering applications. Astrophysical 
occurrences include stellar interiors 111271 |80l l42l l46l and planetary atmospheres and 
interiors Il57ll79l 17^11231 . while terrestrial occurrences include the Earth’s outer core 
BTl [TSl 1^ ri4l . mantle 111 171 . oceans II129I l9Tl . and atmosphere ll^ . Many of these 
systems have internal sources or sinks of buoyancy, including the Earth’s mantle, the 
cores of large main-sequence stars, radiating atmospheres, and nearly any engineered 
system where chemical or nuclear reactions take place in a fluid environment. Among 
such engineering applications, particular attention has been paid to nuclear accident 
scenarios in which exothermic nuclear reactions drive convection in molten material 


sisiisi. 

We speak here in terms of thermal convection, where the body forces are gravita¬ 
tional and depend on the fluid’s density, which, in turn, depends on its temperature. 
Other types of convection are not discussed but are often governed by similar dynam¬ 
ics. In compositional convection, for instance, chemical concentration takes the place 
of temperature. In electroconvection (e.g. ITTIIISIII I. electric charge takes the place of 
temperature, and electrical potential takes the place of gravitational potential. 

Convection can be indefinitely sustained in each configuration studied here, and we 
focus on the time-averaged properties of sustained convection, especially heat trans¬ 
port. Transient phenomena are not addressed. The minimum requirement for ordinary 
fluid to convect is that warmer (less dense) fluid lie below cooler (more dense) fluid, 
and that this adverse temperature gradient be sufficiently destabilizing to overcome 
the viscous forces that damp fluid motion. Eor the convection to be indefinitely sus¬ 
tained against viscous dissipation, there is an additional requirement; an inexhaustible 
source of energy that drives the system away from equilibrium by endlessly adding 
heat somewhere other than the top of the domain or removing heat somewhere other 
than the bottom. This can be accomplished through the thermal boundary conditions, 
as when a pot of water is boiled on a stove, or it can be accomplished through internal 
heat sources or sinks, as when radioactive decay heats the Earth’s mantle. 

The present chapter lays out the basic features of six convective configurations— 
three that are driven solely by the boundary conditions, and three that are driven by 
internal heating. The configurations are defined in M. 1 and the Boussinesq equations 
that are used to model them are given in n.2 followed in §1.3| by our chosen nondi- 
mensionalization. The basic commonalities and differences of the six configurations 
are then summarized: j ]1.4| discusses static states, |1.5| gives a qualitative look ahead 
to the experimental findings that are surveyed in chapter]^ and |1.6| introduces the 
integral quantities and relations that govern heat transport. 


1.1 Six configurations 

The six configurations we study here share the same basic geometry: a horizontal fluid 
layer of height d. In its horizontal dimensions, the layer can be modeled as infinite, 
periodic, or bounded, though only the last case is realizable in physical experiments. 
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Figure 1.1: Schematics of the six configurations studied in the present work. They are 
distinguished by their thermal boundary conditions and by the presence or absence of 
a constant and uniform internal heat source (H). Gravity acts vertically downward. All 
quantities are dimensional; nondimensionalized versions of these schematics appear in 
tables o and |1.2| In the IH2 configuration, H and F are related such that heat loss 
balances internal heat production (see text). 


In a theoretical investigation, we can allow the convection to have three-dimensional 
(3D) freedom, or we can make it two-dimensional (2D) by imposing uniformity in one 
of the horizontal dimensions. 

The six configurations differ only in their top and bottom thermal boundary con¬ 
ditions and in the presence or absence of volumetric heating. Three of the cases are 
versions of Rayleigh-Benard (RB) convection, wherein the flow is driven solely by 
thermal boundary conditions that cause heat to enter across the bottom boundary and 
exit across the top one. The other three cases are instances of internally heated (IH) 
convection. Here, heat is added only by a constant, uniform source, and at least some of 
it exits across the top boundary. The only thermal boundary conditions we employ are 
fixed-temperature, fixed-flux, or perfectly insulating. Fixed temperatures model per¬ 
fectly conductive boundaries, while fixed heat fluxes model boundaries that conduct 
heat poorly EH. 

The thermal boundary conditions of our three RB configurations are shown in the 
top row of figure [LT] In the case we call RBI, which is the most-studied RB model, 
the boundary temperatures are fixed, with a temperature drop of A from the bottom 
boundary to the top one. (The temperature at the top boundary can be fixed at zero for 
convenience because, under the governing equations of our models, only temperature 
differences affect the dynamics.) In RB2, the heat flux across each boundary is fixed 
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by setting the temperature gradients there to the same value, —F. The RB3 case mixes 
the two previous cases, having a fixed-temperature condition on the top boundary and 
a fixed-flux condition on the bottom one. The configuration where these two boundary 
conditions are swapped need not be considered since it is related to RB3 by symmetry, 
at least with the governing equations of our models. 

The thermal boundary conditions of our three IH configurations are shown in the 
bottom row of figure [TTT] Temperature is produced volumetrically at rate H in all three 
cases, corresponding to heat production at rate H jcp, where Cp is the fluid’s heat ca¬ 
pacity. In IHl, the top and bottom temperatures are fixed at the same value, and the 
internally produced heat escapes across both boundaries. In IH2 and IH3, the bottom 
boundary is perfectly insulating, meaning the vertical temperature gradient must vanish 
there. The top boundary, across which all the internally produced heat escapes, has a 
fixed flux in IH2 and a fixed temperature in IH3. In IH2, the boundary flux must match 
the rate of internal heat production in order for convection to be statistically steady, 
so r is determined by H, as described in §1.3| Configurations where the boundary 
conditions in IH2 or IH3 are reversed, making the tops insulating, do not need to be 
considered; fluid in such configurations would remain static, as follows from the sta¬ 
bility results of 12.2 Notice that the temperature is subject to Dirichlet conditions in 
RBI and IHl, Neumann conditions in RB2 and IH2, and one condition of each type in 
RB3 and IH3. 

For the velocity, we impose either no-slip or free-slip boundary conditions. No-slip 
conditions apply to most laboratory experiments and many engineering applications, 
while free-slip conditions are more appropriate in modeling certain astrophysical, geo¬ 
physical, and plasma physical systems. Although we address both possibilities, some 
results are available only for no-slip boundaries. If side boundaries exist, we assume 
they are perfect thermal insulators. If side boundaries exist, we assume they are perfect 
thermal insulators. 

The six models of figure [TT| can be extended in various ways. For instance, inter¬ 
nal heating can be added to the RB configurations, creating hybrid models driven both 
by the boundary conditions and by internal heating. The thermal boundary conditions 
can also be made more complicated, perhaps to model thermal radiation or moderate 
conductivity. However, such models require at least one more control parameter than 
those of figure [TTT] When modeled by the Boussinesq equations, which are introduced 
in the next section, each of the six configurations we study is governed by only two 
dimensionless control parameters, aside from any parameters used to describe the ge¬ 
ometry. In a sense described at the end of [ 1.3 they are the only models of convective 
layers for which this is true, hence they are a natural starting point. 

In many sections of this SpringerBrief results are presented for all six configura¬ 
tions in figure 0 letting us highlight their similarities and differences. Section |2T3] 
and chapter]^ however, focus mainly on IH convection. This work is not meant to be 
a comprehensive review of RB convection, which has been reviewed several times in 
recent decades lll22ll47l [Tl l90ll28]| . Interest in RB convection goes well beyond heat 
transport, as the system has become a canonical model of nonlinear science, having 
provided early examples of instabilities, bifurcations, pattern formation, and chaos in 
spatially extended systems. IH convection, which has been the subject of numerous 
works but is still much less studied than its RB counterpart, was last reviewed in the 
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1980s Il87ll2^ . Bringing this topic up to date requires not only that we review the more 
recent studies of IH convection, which are relatively few in number, but also that we 
reinterpret older studies in light of our contemporary understanding of RB convection. 

Many of the laboratory experiments reviewed in chapter|^are captured well by one 
of our three IH models, but applications are rife with further complications, including 
compressibility, temperature-dependent material properties, complicated geometries, 
chemical and nuclear reactions, rotation, magnetism, and other thermal boundary con¬ 
ditions. Nonetheless, characterizing heat transport in the relatively simple models we 
consider is already a formidable challenge, and the task is far from complete. 

1.2 Boussinesq equations 

A mathematical model of thermal convection must include equations governing the ve¬ 
locity and temperature helds, along with a constitutive relation between temperature 
and density. The compressible Navier-Stokes equations 01401 describe the velocity 
held accurately in a wide range of physical applications, but they are challenging equa¬ 
tions to study analytically or integrate numerically, and they can be avoided when the 
density held does not deviate too strongly from hydrostatic equilibrium. Pressure- 
driven hows with weak density variations can simply be approximated as having con¬ 
stant densities, yielding the incompressible Navier-Stokes equations. In buoyancy- 
driven hows, however, density variations cannot be totally ignored because they create 
the buoyancy gradients needed to drive motion. The typical compromise is to employ 
the Boussinesq approximation, as we do here. 

The Boussinesq approximation, which was hrst invoked by Oberbeck EH and is 
also called the Oberbeck-Boussinesq approximation, involves two main assumptions. 
First, the Huid’s density, p, is assumed to vary linearly with temperature, T, about some 
hydrostatic reference state denoted by p* and T*. That is, 

piT)=p,[l-aiT-T)], (1.1) 

where a is the linear coefficient of thermal expansion. Second, the density variations 
are assumed to be sufficiently weak that they can be ignored everywhere except in 
the buoyancy force. The huid is sometimes called incompressible since the velocity 
held is divergence-free, although compressibility does manifest in the buoyancy vari¬ 
ations. Numerous justihcations have been put forth for replacing the fully compress¬ 
ible Navier-Stokes equations with the simpler Boussinesq equations, typically invoking 
some combination of asymptotic expansions and ad hoc assumptions. See Spiegel and 
Veronis 112611 for one physical justihcation and Rajagopal et al. 010611 for a discussion 
of various other justihcations. The precise assumptions invoked vary, but in all versions 
there is a sense in which gradients of the Huid’s properties should not be too steep. If 
the Boussinesq approximation is used in modeling a physical system, the assumptions 
under which the approximation holds should be checked if possible, either by physical 
measurement or by numerical simulation of compressible equations. 

With constant gravitational acceleration g acting in the —z direction, applying 
the Boussinesq approximation to the compressible Navier-Stokes equations yields the 
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Boussinesq equations lll09lfT8l . 


Vu = 0 

5,u + u- Vu = --^Vp + vV^u + gaTz 
dfT+u-VT = kV^T+H, 


( 1 . 2 ) 

(1.3) 

(1.4) 


where u = (m, v, w) is the fluid’s velocity vector, p its pressure, v its kinematic viscosity, 
and K its thermal diffusivity. The temperature source term, H, is absent from RB 


convection but drives the convection in our IH models. The pressure term in (1.31 has 
absorbed hydrostatic terms of the buoyancy force coming from 


1.3 Nondimensionalization 

To nondimensionalize the Boussinesq equations, we scale distance by the layer height, 
d, time by the characteristic timescale of thermal diffusion, d^/K, and pressure by 
Ptd^fK. We scale temperature by a dimensional quantity. A, that is defined differently 
in various configurations. In the RBI case, A is the prescribed temperature difference 
between the boundaries. In the other cases. 


A;= 


dr RB2, RB3, IH2 
^ IHl, IH2, IH3. 


(1.5) 


Nondimensionalized by these A, the temperature difference between the boundaries is 
unity in RBI; the fixed temperature fluxes are unity in RB2, RB3, and IH2; and the 
volumetric heating rate is unity in the IH cases. Both definitions in ( |1.5| l apply to IH2 
because we add the consistency condition F = dH /k in that case to ensure that heat 
production balances heat loss. 


The Boussinesq equations (1.2 1 -( 1.4 1 in dimensionless form are 


V.u = 0 

(9rU + u. Vu = —VpVPrV^n + PrRTz 
5,7’ + u.V7’ = + 


( 1 . 6 ) 

(1.7) 

( 1 . 8 ) 


where the symbols u, T, p, x, and t henceforth represent dimensionless quantities. The 
vertical extent is 0 < z < 1, and the temperature source term is 


Q 


0 RB 
1 IH. 


(1.9) 


The dimensionless control parameters are the Rayleigh number, R, and Prandtl number, 
Pr, defined by 


R := 


gad^A 


KV 


( 1 . 10 ) 


Pr = 


(1.11) 
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The definition of R differs between cases when the definition ( |1.5| l of A differs. 

The Rayleigh number may be thought of as the ratio of inertial forces to viscous 
forces. When it is large, the fluid is strongly driven by differential buoyancy forces. We 
regard R as the primary control parameter since raising it typically makes the flow more 
complex. For R to indeed be a control parameter, we needed to define the dimensional 
temperature scale. A, in terms of quantities that are known a priori: the boundary con¬ 
ditions or heating rate. However, it is sometimes useful to define a different Rayleigh 
number using a temperature scale that is determined dynamically by the flow. This sort 
of Rayleigh number cannot serve as a control parameter but can be a useful diagnostic 
quantity. Thus, we will sometimes distinguish between th e cont rol Rayleigh number. 


R, and diagnostic Rayleigh numbers, Ra or Ra, defined in 11.6.5 


The Prandtl number is the rate at which the fluid diffuses momentum, relative to 
the rate at which it diffuses heat. Unlike the Rayleigh number, it is a material property 
of the fluid and does not depend on the geometry or boundary conditions. The Prandtl 
number is large in fluids that damp motion strongly and conduct heat poorly. In the 
Earth’s mantle, for instance, Pr is effectively infinite. The Prandtl number is small in 
fluids that damp motion weakly and conduct heat well, such as liquid metals and stellar 
plasmas. Air and water are intermediate examples, having Prandtl number close to 0.7 
and 7, respectively, under atmospheric conditions. 

In all six configurations of figure o the dynamics depend on only two control 
parameters, Pr and R, aside from any parameters describing the geometry, such as as¬ 
pect ratios. This is the minimum number of parameters we can hope for in the study 
of convection, except in those special cases where Pr can be eliminated because it is 
effectively infinite. Additional parameters would be needed if the thermal boundary 
conditions were more complicated 1112511^ . the internal heating law were more com¬ 
plicated ||45111341 . or the boundary conditions and internal heating each introduced 
their own temperature scales llTTl [T9l IfiTl |2l 1124111451 ISTl . We restrict ourselves to 
models that require only Pr and R since every additional parameter makes it much 
harder to understand parameter space. In fact, among the ways that uniform heating, 
fixed-temperature boundaries, and fixed-flux boundaries can be combined, our six con¬ 
figurations (and their symmetry-related siblings) seem to be the only ones governed by 
so few parameters. 

1.4 Static states 


Dimensionless schematics of the RB and IH configurations are shown in the first rows 
of tables o and |1.2| respectively. The basic features of the six cases are summa¬ 
rized in the subsequent rows of both tables and are further laid out in the remainder 
of this chapter. The simplest solutions to the governing equations are static, with heat 
transported only by con duction. These are the unique asymptotic states when R is suf¬ 
ficiently small (cf. ^ 2.2 1 , and they solve the Poisson or Laplace equation V^T + Q = 0. 


Since we assume that side boundaries are nonexistent or perfectly insulating, the static 
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^(r)^ i(i-(w7’)) 
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1 - {wT) 


Table 1.1: Summary of the properties of RB convection discussed in the present chap¬ 
ter. All quantities are dimensionless, and the vertical extent is 0 < z < 1. Notation is 
dehned throughout the chapter. Briefly, * denotes an average over horizontal directions 
and time, (*) denotes an average over volume and time, 5(7’) is the mean temperature 
of the fluid relative to that of the top boundary, and J = —d^T + wT is the sum of the 
conductive and convective vertical heat fluxes. 
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Table 1.2: Summary of the properties of IH convection discussed in the present chapter. 
All quantities are dimensionless, and the vertical extent is 0 < z < 1. Notation is defined 
throughout the chapter. Briefly, * denotes an average over horizontal directions and 
time, (*) denotes an average over volume and time, d{T) is the mean temperature of 
the fluid relative to the top boundary, J = —d{T + wT is the sum of conductive and 
convective vertical heat fluxes, and Tmax is the maximum value of ^(z). 
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temperature fields, Tgt, vary only in z: 


1 —z 

Ztiz) = <( iz(l -z) 
. 5 ( 1 -^") 


RBI, RB2, RB3 
IHl 

IH2, IH3. 


( 1 . 12 ) 


These purely conductive profiles are depicted in the second rows of tables [LT] and 1.2 


They are parabolic with internal heating and linear without it. The static profiles in 
RB2 and IH2 are determined only up to additive constants, but these constants do not 
affect the dynamics, so we have fixed them for convenience. 

In each configuration, our dimensional temperature scale. A, is characteristic of 
the static state. This is why the dimensionless T, have no dependence on R. When 


we define diagnostic Rayleigh number in 11.6.5 we will do so by replacing A with 
temperature scales characteristic of the convective states, rather than the static ones. 


1.5 Temperature fields in strong convection 

Whereas the fluid remains static when R is sufficiently small, it convects strongly when 
R is large. Convection typically strengthens monotonically as R is raised, though 
this is not universally true and can depend on how strength is quantified. (The non¬ 
monotonicity of convective transport in 14^ provides a counterexample.) Figure [L2] 
shows instantaneous temperature fields from 2D simulations at large R. The RB1 field 
in figure [L2^ is representative of all three RB cases; hot plumes rise from the bottom 
boundary, cold plumes descend from the top one, and both types of plumes contribute 
to upward heat transport. In the IHl field of figure cold plumes descend from 
the top, but the bottom boundary layer is cold and stably stratified. This bottom layer 
emits no buoyant plumes, so any mixing between it and interior must be driven by 
shear, rather than buoyancy. The IH3 field in figure |1.2{ : is representative of both IH2 
and IH3; cold plumes descend from the top boundary, and there is no thermal boundary 
layer at the bottom. 

The turbulent convection that occurs at large R creates mean vertical temperature 
profiles very different from the static ones. Rough schematics of such profiles are 
shown in the third rows of tables o and |1.2| Although these schematics include no 
secondary details, they illustrate the main differences between the various configura¬ 
tions. In each of the six cases, mixing by strong convection tends to flatten the tempera¬ 
ture profile in the layer’s interior. (We are not aware of a counterexample in 3D, though 
one exists in 2D under conditions that allow very strong winds to develop 14^ 11041 .) 
The roughly isothermal interiors are flanked by one or two thermal boundary layers, 
and these are what distinguish the various cases. 

In all six of our models, temperature is unstably stratified at the top boundary. 
At the bottom boundary, the temperature is unstably stratified in the RB cases, sta¬ 
bly stratified in IHl, and unstratified in IH2 and IH3. These facts are evident in the 
static temperature profiles (cf. tables o and [13 and remain provably true of sus¬ 
tained convection, at least in a time-averaged sense. As convective heat transport rises, 
the mean temperature profiles undergo various changes. In RBI, where the tempera¬ 
ture difference between the boundaries cannot change, the boundary layers steepen as 
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Figure 1.2: Temperature fields from 2D simulations of (a) RBI, (b) IHl, and (c) 1H3. 
Each simulation employed a horizontal period of 6, no-slip boundaries above and be¬ 
low, Pr = 1, and R/Rl = 10^, where Rl is the Rayleigh number at which the static 
becomes linearly unstable (cf. ^2.1 1 . The coolest fluid (blue) has a temperature of zero 
in each case, and the warmest fluid (red) has a temperature of 1, 0.017, and 0.044, 
respectively. 
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the heat flux through the domain rises. In the other two RB cases, where the mean flux 
through the domain cannot change, the temperature difference between the boundaries 
decreases. In the IH cases, the temperature of the fluid, relative to that of the top bound¬ 
ary, drops as the convection strengthens. The produced heat leaves only across the top 
boundary in IH2 and IH3. It leaves across both boundaries in IHl, but the majority 
leaves across the top boundary, hence the top boundary layer is steeper than the bottom 
one. 

Although the turbulent temperature profiles are fairly easy to understand qualita¬ 
tively, it is very difficult in general to anticipate their quantitative features. This would 
be tantamount to accomplishing our main objective of characterizing the bulk heat 
transport. 

1.6 Mean heat fluxes and integral relations 

Heat in a convecting fluid is transported by two mechanisms simultaneously: conduc¬ 
tion and convection. Conduction refers to the diffusion of heat down the temperature 
gradient, while convection refers to the advection of heat by fluid motion. The tem¬ 
perature equation ( |1.8[ ) can be written in the standard form of a conservation law as 
dfT + V ■ J — Q, where J := uT — VT is the total heat current at a point. Evidently, the 
conductive current is —VT in our nondimensionalization, and the convective current 
is uT. The horizontal components of J vanish when averaged over the volume since 
our side boundaries are insulting or nonexistent. Here we focus on the heat current’s 
vertical component, J, where 


(1.13) 

(1.14) 



Much of our effort is devoted to quantifying the relative contributions to vertical heat 
transport made by conduction and convection—that is, by —d{r and wT. 

In our notation, an overbar, as in /, denotes an average over horizontal surfaces and 
infinite time. Angular brackets, as in (/), denote an average over the entire volume and 
infinite time. When the dimensionless domain is bounded horizontally by 0 < .r < L;c 
and 0 < y < Ly, 



(1.15) 


(1.16) 


The above limits can be replaced with liminf or limsup to ensure their existence. For 
simplicity, we assume in our calculations that infinite-time averages commute with ver¬ 
tical averages and that horizontal averages commute with vertical derivatives, though 
these assumption can often be avoided. In the above notation, the mean heat flux across 
a horizontal surface is 


J{z) = -T'{z)+wT{z), 


( 1 . 17 ) 
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where the prime indicates differentiation in z- The mean vertical flux across the entire 
layer is 

{J) = 5T + {wT), (1.18) 

where the mean conductive flux, 

dT:=TB-TT, (1.19) 


is the difference b etwee n the mean bottom temperature, Tb, and mean top tempe rature, 
TT- Expressions (|l.l7| and ( 1.18 1 are simply averages of the defini tion ( 1.141 for 7; 


conflguration-speciflc constraints on J{z) and (7) are given in [ 1.6.2 


1.6.1 Heat balances 


Conservation of heat energy is expressed in the various cases by the heat balances 

RBI, RB2, RB3 (1.20) 

IHl (1.21) 


t' — t' 

It — ^ R 


-t' 


- t '^1 


-T't = 1 


IH2, IH3. 


( 1 . 22 ) 


These balances, which are shown also in the fourth rows of tables [TT| and |1.2| are de¬ 
rived by averaging the temperature equation (1.8 i over volume and time. Time deriva¬ 
tives vanish from such averages because the instantaneous volume averages of u and 


T are bounded uniformly in time (cf. (2.31. In the RB cases, the balance reflects the 
fact that mean heat flux into the layer at the bottom (—Tg) must equal the mean flux 
out of the layer at the top (—Tj). In IHl, the rate of internal heat production (unity) is 
balanced by the combined outward fluxes of heat across the top boundary {—Tj) and 
the bottom one (T^). In IH2 and IH3, where the bottom boundary is insulating, the 
internal production is balanced entirely by the outward flux across the top boundary 
(-T't). 


1.6.2 Constraints on net heat fluxes 


Little can be said a priori about the variation with height of the mean heat flux compo¬ 
nents, — T (z) and wT{z), but we can derive constraints on their sum, 7(z). Averaging 
the temperature equation (|1.8|l horizontally, vertically from 0 to z, and temporally gives 


J{z) := —T'{z) + wT (z) 


-T'g RBI 

1 RB2, RB3 

-t'b+z IHl 

z IH2, IH3. 


(1.23) 


The net vertical flux is the same at every height in the RB cases and increases linearly 
with height in the IH cases. In the four cases where a boundary flux is fixed, 7(z) is 
known exactly. In RBI and IHl, 7(z) is determined only up to the mean heat flux at 
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the bottom boundary, —T'q. Below we give some alternate expressions for J{z) in these 
two cases, but they all involve quantities that, like —Tb, are not known a priori. 

The volume-averaged heat flux, (7), is also constrained. In th e four cases where 
J{z) is known exactly, (7) is found by vertically integrating (1.23 i. In RBI and IHl, 
the fixed-temperature boundary conditions ensure that the mean conductive fluxes are 
5T = \ and 5T = 0, respectively. The results are 


(7) :=5T + {wT) = 


l + {wT) RBI 

1 RB2, RB3 

(wT) IHl 

i IH2, IH3. 


(1.24) 


In all six cases, we would like to know how the control parameters affect the convective 
flux, {wT). In RBI and IHl, where 5T is fixed, this is equivalent to knowing (7). In 
the other four cases, where (7) is fixed, it is equ ivalen t to knowing dT. 

In the RBI and IHl cases, the expressions (1.231 for J{z) can be rewritten by re¬ 


placing Tg with other int egral q uantities. Relations between T^ and Tj are provided 
by the heat balances of (1.6.1 while relations between t'b and {wT) are found by 


equati ng the (7) expressions (1.24 1 with the vertical integrals of the 7(z) expressions 
(1.23 1 . The alternate expressions for 7(z) found in this way are 


J(z) = 


—t' 

^ B 

-t'b+z 


— t' 

1 j 


= -Tr-il-z) 


= l + {wT) 

= D + * 


vT) 


RBI 

IHl. 


(1.25) 


The constraints on 7(z) and (7), expressed in terms of volume integrals, are summa¬ 
rized in the fifth and six rows of tables ll.ll and [L2l 

In the IHl configuration, there is yet another useful way to interpret (wT): in terms 
of the fractions of internally produced heat that flow outward across the top and bottom 
boundaries. Expressions for these fractions, which we call and follow from 
relations p.21| ) and ( |1.25| ), 

=-f y = 1-P (wT) (1.26) 

t'b=\-{wT). (1.27) 

The top and bottom fractions are both 1/2 in the static state, but convective transport 
breaks this symmetry. 


1.6.3 {wT) and 8{T) 

Essential information about heat transport is captured by the volume integrals {wT) 
and 5(7’), where 

5{T) := {T-Tt) (1.28) 

is the mean fluid temperature, relative to that of the top boundary. The above definition 
is needed only for RB2 and IH2, where the top temperature is not fixed. In the other 
configurations, 5(7’) = (T) since we have set Tt = 0. Neither (wT) nor 5(7’) is known 
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a priori when the fluid is flowing. Instead, the quantities must be studied by physical 
and computational experiments, and they can sometimes be bounded analytically. All 
bounds and many experimental findings that we review in the following chapters can 
be stated in terms of (wT) or 5(7’). In the literature, however, results on (wT) are often 
stated differently but equivalently in terms of other quantities, including 5T, Tg, Tj. 


and the Nusselt number N defined in ^ 1.6.4 below. 


1.6.3.1 Uniform bounds 


The seventh and eighth rows of tables o and |1.2| give uniform bounds on (wT) and 
5(7’)—that is, bounds that are independent of R and Pr. The bounds are derived in this 
chapter’s appendix. The lower bounds on (wT) are tight in all six configurations. The 
upper bounds are thought to be tight among uniform bounds, except in the IH2 case, 
where we suspect a uniform upper bound of 1/2. In the IH cases, the upper bounds on 
5(r) are tight, and the lower bounds are thought to be. In the RB cases, on the other 
hand, it is not clear whether any of the bounds on 5(7’) are tight. 

The mean convective flux, {wT), saturates its lower bound of zero in each con¬ 
figuration only in the static state. Physically, this is because {wT) is proportional to 
the work exerted by buoyancy, and when motion persists this work must be positive to 
balance viscous dissipation. Mathematically, the positivity of (wT) in sustained con¬ 
vection follows from relation (2.37 1 in the next chapter. The upper bounds on (wT) 


correspond to limits in which convective transport is infinitely stronger than conduc¬ 
tive transport. In the RBI case, where 57’ = 1, this limit is approached when (wT) 
grows without bound. In the four cases where the total heat flux, (7), is fixed, this limit 
is approached as (wT) —?> (J) and 57’ 0. The IHl case is different in that 5T — 0, 

so (wT) is solely responsible for the mean vertical flux. However, if one thinks of the 
outward transport of heat across both boundaries, rather than upward transport, then 
the upper limit (wT) -^1/2 indeed means that convection fully takes over from con¬ 
duction. The corresponding limits of the top and bottom flux fractions ( |1.26| l-( [TT27] i 
are > 1 and —>■ 0. 

The mean temperature relative to the that of the top boundary, 5(7’), is bounded 
above and below in all six cases, but the IH bounds differ in character from the RBI 
bounds. The RBI bounds given in table o ensure that the mean temperature of the 
fluid lies between those of the top and bottom boundaries. The same may be true of 
5(7’) in RB2 and RB3, but the bounds derived in this chapter’s appendix are too crude 
to show it. In RB convection, the mean fluid temperature is exactly halfway between 
the boundary temperatures in the static state. The same is often true when the fluid is 
flowing, at least with symmetric boundary conditions, but it seems no rigorous state¬ 
ments have been proven that reflect this observation. In the IH cases, 5(7’) saturates 
the upper bounds of table |1.2| only in the static states and is strictly smaller when the 
fluid is flowing. Its lower bound of zero, much like the uniform upper bounds on (wT), 
corresponds to convection being infinitely stronger than conduction. The 7?-dependent 


bounds of) 2.3 show that 5(7’) could approach zero only as 7? —?> oo. 

In IH convection, where R is proportional to the dimensional heating rate, 77, it 
might seem counterintuitive that raising R tends to decrease 5(7’). However, the di¬ 
mensional mean temperature, 5(7’) A, indeed rises with 77. The dimensionless quantity 
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RBI 
RB2, RB3 
IHl, IH2, IH3 


R-dependent R-dependent 
bound on {wT) bound on 5{T) 
(wT) < none 

{wT) < 1 -cR^^/ 3 none 

none 5(7’) >cR^*/^ 


Table 1.3: Proven bounds on {wT) and 5(7’) that hold at large R. The constant c differs 


between cases. Details and references are given in ) 2.3 


5{T) falls as convection strengthens because it has essentially been normalized by its 
static value. 


1.6.3.2 R-dependent bounds 


Since many of the uniform bounds on {wT) and 5(7’) given in tables o and|1.2| 


tight, improving them requires finding bounds that depend explicitly on R or Pr. Some 
R-dependent bounds have been proven for the configurations we are considering, as 
summarized in table 1.3 Bounds that vary with Pr have been proven recently for RBI 
ll29l but not yet for other cases, although some bounds have been proven for the infinite- 
Pr limit that are tighter than the corresponding uniform-in-Rr results, as discussed in 


In RB convection, the R-dependent bounds that have been proven are all upper 
bounds on {wT). They approach the uniform upper bounds as R —oo but are tighter at 
all finite R. The uniform lower bounds of 0 < {wT) cannot be improved upon, if they 
are to hold for all solutions, since they is saturated by the static states. These states 
are unstable at large R, however, and {wT) typically grows with R in experiments and 
simulations. There might exist better lower bounds that hold only for attracting states, 
rather than all solutions, but we lack the mathematical machinery to prove them. 

In IH convection, the R-dependent bounds that have been proven are all lower 
bounds on 5{T). They approach the uniform lower bounds of zero as R —oo but are 
tighter at all finite R. No R-dependent upper bounds on {wT) have been proven in IH 
convection, but we argue in |1.6.3.4| that such proofs should be possible. On the other 
hand, the uniform upper bounds on 5(7’) and lower bounds on {wT) are saturated by 
the static states, so any efforts to improve them run into the same obstacle as efforts to 
improve the lower bounds on {wT) in RB convection. 


1.6.3.3 Empirical approximate relations 

In addition to the exact integral relations and bounds discussed above, experiments and 
simulations suggest some approximate relations for 5(7’) at large R. These relations are 
summarized in the ninth rows of tables [TT| and |1.2| Most can be expressed as relations 
between 5(7’) and {wT), at least approximately, but the underlying assumption in IH 
convection differs from that in RB convection. 
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In the RB cases, the underlying assumption is that the top and bottom boundary 
layers are roughly symmetric. This implies that the mean fluid temperature is about 
halfway between the top and bottom temperatures, as in the schematics of turbulent 
temperature profiles in table|l.l| That is. 


5{T) - \5T^ 


{l-{wT}) 


RBI 

RB2, RB3 


(1.29) 


for large R. These relations are not expected to hold exactly when the top and bottom 
boundary conditions differ, but they could nonetheless be approached as R —> oo. It 
might be possible to prove precise versions of the above statements, such as upper and 
lower bounds on 5{T) that converge to 5T as R^oo, but we are not aware of any such 
results. 

In the IH cases, the underlying assumption is that mean temperature profile, T{z), 
at large R is roughly isothermal outside of thin thermal boundary layers, as in the 
schematics of turbulent temperature profiles in table 1.2 Experimental support of this 
assumption is presented in chapter]^ Approximate isothermally in the IH cases implies 
that 


5{T)^ 


5r = i- 


•wT) 


IHl 

IH2, IH3 


(1.30) 


for large R, where Tmax is the maximum mean temperature that T (z) attains in the layer. 
In IHl, the assumption of an isothermal interior does not give a relation between 5{T) 
and {wT), nor is any simple relation suggested by experiments. 


1.6.3.4 Conjectured upper bounds on {wT) in IH convection 

In IH2 and IH3, the empirical observation that 5{T) ^ 5T at large R suggests that the 
two quantities might obey similar bounds. Since bounds of the form 5{T) > cR^^/^ 
have been proven, it seems likely that bounds of the form 5T > cR^ could be proven 
also. The latter can be alternately stated as upper bounds on {wT): 

Conjecture 1 In the IH2 and IH3 configurations, there exists a constant c > 0 such 
that for all Pr and sufficiently large R, 

(wT) < 

In IHl, experiments suggest that the growth of (wT) with R is similarly bounded 
above (cf. chapter]^. However, since no empirical relation between 5(7’) and (wT) is 
apparent, the proven lower bound on 5(7’) does not suggest a form for an upper bound 
on {wT). We speculate that the upper bound should approach the uniform bound of 
1 /2 algebraically as R —> oo, but we cannot anticipate the algebraic power; 

Conjecture 2 In the IHl configuration, there exist constants c > 0 and a > 0 such 
that for all Pr and sufficiently large R, 

(wT) < 5 -cR^“. 
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1.6.4 Nusselt numbers 

The relative strengths of convective and conductive heat transport are often quantified 
using dimensionless Nusselt numbers. In RB convection, the typically used definitions 
of Nusselt numbers can all be expressed in terms of (wT). In IH convection, dimen¬ 
sionless quantities resembling the RB Nusselt numbers can be defined in various ways 
by invoking {wT), 5{T), or Tmax- Here we consider two ways of defining Nusselt- 
number-like quantities, N and N. The quantity N is determined by T,„ax in IHl and 
by {wT) in the other five cases, while the quantity N is determined in the IH cases by 

5(r). 


1.6.4.1 The Nusselt number N 


The definition of N that we choose is one that has helped reveal parallels between 
various RB configurations when used in concert with the quantity Ra defined in the 
next subsection 199111471 17^ 11531 . In every case other than IHl, our definition of N 
can be expressed as the ratio of mean total heat flux to mean conductive heat flux, 
where both quantities are averages over volume and time in the developed flow. 


N = 


(J) 5T + {wT) 


{Jcond) 


5T 


RB1,RB2, RB3, IH2, IH3. 


(1.31) 


The above definition would fail for IHl because its deno minat or would be zero. Ap¬ 
plying the various constraints on ST and (wT) (cf tables and|l.2| to expression 
( 1.31|l and adding an ad hoc definition for the IHl case, we obtain 


N\= 


' 1 + (wT) 

1 _ 1 

Wf ~ 1 - (wT) 

1 


^T,fui. 

1 


1 

l-2(w7’) 


RBI 

RB2, RB3 
IHl 

IH2, IH3. 


(1.32) 


The rationale for our definition of N in the IHl case, the only case where heat 
flows outward across both boundaries, is that we are considering outward heat fluxes 
instead of upward fluxes. The mean total outward flux is unity since it must balance 
heat production. To determine the mean outward conduction, we imagine dividing the 
layer at the height z* where the temperature profile T (z) assumes its maximum value 
of Tmax- The upward conduction above z* is proportional to Tmax, as is the downward 
conduction below z* ■ Thus, the ratio of total outward transport to conductive outward 
transport is inversely proportional to Tmax- The 1 /8 factor makes N unity in the static 
state. Although the analogy between A in IHl and in the other five cases is not perfect, 
the experiments reviewed in chapter[^reveal similarities in N between all cases. In the 
IH3 case, the quantity we call N has been considered under various names, perhaps 
first by Thirlby 111391 . In the IHl case, our definition has apparently not been used, but 
many authors have considered T,nax, as well as the so-called top and bottom Nusselt 


numbers discussed in ^ 3.2.2 
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In all cases except IHl, our definitions (1.321 of N can be expressed in terms of 


(wT) alone. One might wonder whether N in the IHl case would be better defined as 
inversely proportional to 1 —2{wT) instead of to Tmax- This would be superficially 
identical to the IH2 and IH3 definitions of N, provided the la tter are expressed in terms 
of {wT). However, the IHl experiments discussed in ^3.2 confirm that T,„ax behaves 


very much like an inverse Nusselt number, while the quantity 1—2 {wT) does not. 


As described at the end of ^ 1.6.2 above, {wT) in IHl convection instead conveys the 


asymmetry between upward and downward heat fluxes. 


1.6.4.2 The Nusselt number N 

Since 5(7’) is physically important in IH convection, it is natural to define a Nusselt- 
number-like quantity that is exactly related to 5{T), rather than to (wT) or T„,ax- It 
works well to simply define N as inversely proportional to 5(7’), 


1 


N:= 


125(7’) 

1 

,35(7’) 


IHl 


IH2, IH3. 


(1.33) 


In the IH2 and IH3 cases, this definition could be interpreted as 


N = 


cond ) 


IH2, IH3, 


(1.34) 


which is like the expres^on ( |1.31| l for N with averages weighted proportionally to 
height. We do not define N for RB convection, although the mean temperature in those 
cases merits attention also, as discussed in §1.6.3| 


1.6.4.3 Basic properties 

In almost all cases it has been proven that N > I and N > I, with equality holding only 
in the static states. These facts follow from the uniform bounds on {wT) and 5(7’) 


discussed in ) 1.6.3.1 It remains to be proven that IV > 1 in the IHl case, which would 
be true if Tmax could not exceed its static value of 1/8. In turbulent convection, it is 
typically expected that N ^ and N ^ In the IH cases, this is tantamount 

to Tmax —> 0 or 5(7’) —> 0. Such limiting behavior has not been proven but is supported 
by the ex perim ental results described in chapter]^ 


Table 


1.4 


summarizes the /^-dependent bounds that are known for N and N. These 

In 


are simply restatements of the bounds on {wT) and 5{T) given above in table 1.3 
RB convection, the upper bounds on N are equivalent to upper bounds on {wT). In 
IH convection, the upper bounds on N are equivalent to lower bounds on 5(7’). Upper 
bounds on N have not yet been proven for IH convection. In IH2 and IH3, bounds of the 
for m IV <c R^I^ would follow from the upper bounds on {wT) that we have conjectured 
in ) 1.6.3.4 A bound of the same form for IHl would require showing that T^ax decays 
no faster than 7?^'/^. The quantity Tmax seems harder to access mathematically than 
volume averages like 5(7’) and {wT), which arise naturally in integral relations. 
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RBI 
RB2, RB3 
IHl, IH2, IH3 


R-dependent R-dependent 
bound on {wT) bound on 5{T) 
N < cR^/^ none 

N <cR^/^ none 

none N <cR^I^ 


Table 1.4; Proven R-dependent bounds on N and N that hold at large R. The constants c 
differ between cases. These are re-expressions of the bounds on {wT) and 5(7’) shown 
in table ll.3l 


1.6.5 Diagnostic Rayleigh numbers 

The primary purpose of defining N and N as we have is to identify similarities between 


the various configurations. The bounds in table 1.4 suggest that we have almost suc¬ 


ceeded, but the RBI exponent is 1/2, while the others are 1/3. To bring the various 
cases completely into alignment, we must speak of the dependence of N and N on di¬ 
agnostic Rayleigh numbers, Ra and Ra, instead of on the control Rayleigh number, R. 
These diagnostic parameters can be written simply as 


Ra : 


R RBI 

R/N RB2,RB3,1H1,IH2, IH3 
lla-.^R/N IH1,1H2,1H3. 


(1.35) 

(1.36) 


In terms of these variables, the RB bounds in table 

- 1/2 


1.4 


all take the form N < cRa^!'^, 
and the IH bounds take the form N < cRa ' . Moreover, the analogies brought out by 
considering N and Ra (or N and Ra) are not limited to bounds; similarities emerge also 
in experimental data 0147117^ and heuristic scaling arguments lISTIl . 

The different definitions of R, Ra, and Ra can be viewed as differences in the tem¬ 
perature scale used to defin e a Rayleigh number. The dimensional temperature scales 
A used to define R in ) 1.3 are characteristic of the static states, whereas Ra and Ra 
effectively replace A with temperature scales of the flowing fluid, Aru and A^. In IHl, 
the temperature scale of Ra is the maximum horizontally averaged temperature in the 
flowing fluid. In the other five cases it is the mean temperature difference between the 
boundaries in the flowing fluid. That is. 


^Ra = 


A 

RBI 

5TA 

RB2, RB3 

^Tmax^ 

IHl 

25TA 

IH2, IH3. 


(1.37) 


Replacing A with Arq in the definition (l.lOi of R leads to the above definition (1.35 i 
of Ra. In the IH cases, the temperature scale of Ra is the volume-averaged temperature 
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of the flowing fluid, 


-‘fia 


125(7’)A IHl 
3d{T)A IH2,IH3. 


(1.38) 


Replacing A with A^ in the definition (l.lOi of R leads to the above definition (1.36 1 


of Ra. 


Appendix 

In this appendix we prove various bounds on (wT) and 5{T) that are uniform in the 
parameters, R and Pr. Most of these results are standard, but it is difficult to trace their 
origins, and we do not try. 

Extremum principles 

In each configuration with T = 0 on a boundary, there holds a minimum principle 
giving pointwise, instantaneous lower bounds on T{x,t). For simplicity we assume that 
solutions to the Boussinesq equations exist and remain smooth. If T (x,f) ever achieves 
a local minimum on the interior, then at that point u • Vr = 0 and V^T > 0, and so 
dfT >0. Thus, if the interior is initially warmer than the fixed boundary temperature of 
zero, it remains warmer for all time. Even if part of the interior is initially cooler than 
the boundary, it will be warmer than the boundary at large times. In the RBI case, an 
analogous maximum principle holds relative to the warmer boundary, on which T =\. 
For all X on the interior and sufficiently large f, 

7’(x,f)>0 RB1,RB3, IH1,IH3 (1.39) 

r(x,f) < 1 RBI. (1.40) 

In the RB2 and IH2 configurations, where fixed-flux conditions are imposed on both 
boundaries, neither maximum nor minimum principles hold pointwise. 


Mean convective transport 


Uniform bounds on the mean convective flux, {wT), in our RB and IH configurations 
are summarized in tables o and |1.2| Many of these bounds follow from the power 
integrals (|2.37|i-(|ZM|, which are stated below for convenience. 


(|Vu|2)=R(wr) 

ri + (w7’) RBI 

(ivrl^) = i 5r = l-(w7’) RB2, RB3 

[5{T) IH1,IH2, IH3. 


In all six configurations, the u power integral implies (wT) > 0. Since (|Vu P>> Oif 
convection persists, (wT) saturates its lower bound of zero if and only if the system 
approaches the static state as f —?► 0 °. 
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The uniform upper bounds on (wT) given in tables [0] and |1.2| follow in most cases 
from lower bounds on 5T. In RBI there is no upper bound on (wT). We get 5T > 0 
from the T power integral in RB2 and RB3 and from the minimum principle in IH3. 
This lower bound on ST gives {wT) < 1 in RB2 and RB3, where {wT) + 5T — and 
it gives {wT) < 1/2 in IH3, where {wT) -\-5T = 1/2. These upper bounds on {wT) are 
probably approached by certain solutions, including the turbulent attractors, as R —> oo. 
If so, they are tight among uniform bounds. In IHl, the upper bound {wT) <1/2 
follows from the minimum principle since 1/2 — {wT) = Tg >0. 

It is likely that 5T > 0 in IH2 also, but we settle for the cruder estimate 5T > 
— 1 /V^, derived as follows. 

|5t| =1(5,7')! 

<(l^zr|) 

<{\vt\Y^ 

-A- 

The third line above follows from the Cauchy-Schwarz inequality, the fifth line follows 
from the T power integral, and the last line follows from the bound 5(7’) <1/3 that is 
proven below. The inequality is in fact strict since 5(7’) <1/3, except in the static state. 
Since (wT) + 57’ = 1/2 in IH2, the lower bound on 57’ gives (wT) < 1/2+1 /V^. 

In the IHl configuration, we have also discussed T^ax, the maximum value that 
T{z) attains. For this quantity, the lower bound T,„ax > 0 follows from the minimum 
principle, and the upper bound 7’mm- < 1 / follows from a calculation very similar 
to the one in the previous paragraph. This upper bound is likely not tight; it might well 
be that T^ax never exceeds its static value of 1/8. 


Mean temperature 

Uniform bounds on the mean fluid temperature relative to that of the top boundary, 
5(r), are summarized in tables[TT|and 1.2 In the IH cases, the lower bounds d{T) > 0 


follow from the T power integral. The upper bounds are proven by integrating 


against the T equation (|1.8[), and using ( |1.26| l in the IHl case, to find 


5(7’) = 


IHl 

\ - (zwT) IH2, IH3. 


(1.41) 


Incompressibility gives w = 0 and thus {wT) = {w9) and {zwT) = (zwO), where 9 is 
the devia tion of T from its static profile. Integrating the temperature fluctuation equa¬ 
tion (2.3 I against 9 gives ((z — j) w0) = (| V0|^) > 0 in IHl and {zw9) = (|V0|^) > 0 
in Illland IH3. Therefore, 5(7’) < 1/12 in IHl, and 5(7’) < 1/3 in IH2 and IH3. 

In the RB cases, none of the uniform bounds on 5(7’) are likely to be tight. The 
extremum principles give 0 < 5(7’) < 1 in RBI and 0 < 5(7’) in RB3. The upper 
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bound for RB3 and the upper and lower bounds for RB2 follow from the inequality 
|5(r) I < 1/ that is derived below. 

|5(7')|=|(z5,r)| 

<(k^zr|) 

The first line of the derivation follows from integration by parts, the third line follows 
from the Cauchy-Schwarz inequality, the fifth line follows from the T power integral, 
and the last line follows from the bound 5T = 1 — (wT) < 1 that is proven above. 



Chapter 2 
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The preceding chapter has defined six convective configurations and described, for 
each case, how various integral quantities characterize bulk heat transport. The quan¬ 
tities of central importance include the mean vertical transport by convection, (wT), in 
RB and IH convection and the mean temperature of the fluid relative to that of the top 
boundary, 5(7’), in IH convection. We would like to predict the values that these quanti¬ 
ties assume for various Rayleigh numbers, Prandtl numbers, and confining geometries. 
This is equivalent to predicting the parameter-dependence of the Nusselt numbers we 
have defined; N is determined by T,„ax in IHl and by {wT) in the other five cases, and 
N is determined by 5{T) in the three IH cases. The task before us is very difficult in 
general and would require a greatly improved understanding of fluid turbulence, so we 
are limited to partial results. The present chapter presents facts about heat transport 
that can be determined purely by mathematical analysis of the Boussinesq equations, 
while the next chapter addresses simulations and laboratory experiments. 

There are two main ways of studying heat transport analytically; examining simple 
particular solutions that can be written down exactly or asymptotically, and deriving 
bounds on integral quantities that apply to all solutions. The first method yields much 
stronger results but is useful only at small Rayleigh numbers, where the system either 
remains static or assumes a simple flow. Here we consider only the static states. Heat 
transport in static states is purely conductive and easy to understand, so the main task 
is to determine the parameters at which such states are stable. To this end, we can find 
a Rayleigh number, R^, above which a static state is linearly unstable, and a Rayleigh 
number. Re, below which we can prove that it is the unique globally stable state. These 
results work together with the second method of analysis—bounding N or N above 
by functions or Nb{R )—to constrain the dependence of Nusselt numbers on R. 

Sections 


2.U2.3 


I outline the calculations and values of Re, Re, and Nt,{R), respectively, 
for the various configurations. 


The schematic of figure 2.1 shows how Re, Re, and Nb{R) combine to give some 


knowledge of N in IH convection. In the lowest-/? regime, where R < Re, we know 
that N = 1. This is because the system asymptotically approaches the static state, 
so its Nusselt number, which we have defined as an infinite-time limit, must be that 
of the static state. In the subcritical regime, where Re < R < Re, the static state is 
linearly stable, but sustained flow might also be possible, so all we can say is that 
1 < N < Nb{R) in this regime. In the larger-R regime where Re < R, the static state 
is linearly unstable, so any physically realizable state must have sustained flow and 
thus a Nusselt number strictly greater than unity. That is, 1 < iV < Nb{R) for attracting 
states, although N = I remains possible if unstable states are included. A schematic 
like figure [ZT| for the other Nusselt number, N, in IH convection would lack the upper 
bound since /^-dependent bounds on {wT) have not yet been proven (cf. ) 1.6. 3.4 1 . 

Figure [2T| represents a scenario where Re is strictly smaller than Re- In RB convec¬ 
tion there is no subcritical regime since Re = Re- This allows for asymptotic solutions 
in the weakly supercritical regime, giving more precise expressions for Nusselt num¬ 
bers there. In IH convection, on the other hand, subcritical convection is not ruled out 
because Re < Re- Stronger methods of analysis may be able to prove stability thresh¬ 
olds larger than Re, but not necessarily as large as Re- Subcritical convection is indeed 
possible in IHl 01431 and IH3 0118111441 fT^ . In IH2, the possibility of subcritical 
convection remains open. 
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R. 


Figure 2.1; Schematic of what this chapter’s analytical results say about the depen¬ 
dence of N on R in IH convection. The numerical values of Re, Rl, and Nb{R) vary 
between conhgurations. In RB convection, the analogous diagram for lacks a middle 
region since Re = Rl- 



RBI 

RB2 

RB3 

IHl 

IH2 

IH3 

Re 

1707.76 

720 

1295.78 

26926.6 

1429.86 

2737.16 

Rl 

1707.76 

720 

1295.78 

37 325.2 

1440 

2772.27 

Nb 

Nb 

0.027 

0.28Ra'/2 

0.28 

none 

0.025^'^^ 

none 

0.13 

none 

0.094 


Table 2.1: For no-slip boundary conditions: a Rayleigh number above which the static 
state is linearly unstable {Rl), a Rayleigh number below which the static state is glob¬ 
ally stable {Re), and upper bounds {Nt and Nj,) on the Nusselt numbers N and N that 
are valid for asymptotically large R. References for these values are given throughout 
the chapter. In the IH cases, R-dependent bounds on N have not been proven. In the 
RB cases, N is not dehned. 


For no-slip c ond itions on the velocity, hgure[ZT]is made concrete in each conhg- 


uration by table 2.1 which gives values for Re, Re, Nb{R), and A4(R). The bounds 


have been simplihed by assuming that R is asymptotically large, and they are stated in 
terms of the diagnostic Rayleigh numbers, Ra and Ra, that we have dehned in §sec: 
diagnostic Rayleigh numbers. (Recall that Ra equals R in RBI but equals R/N in the 
other hve cases, and that Ra = R/N in IH convection.) The similarities between the 
various bounds are evident. The present chapter explains how the values in table |2.1| 
are calculated and gives values for some other boundary conditions on the velocity. 

The linear and nonlinear stability analyses that we apply to the static state can be 
applied to other particular solutions as well. Such analyses must be carried out asymp¬ 
totically or numerically, however, since none of the hnite-amplitude particular solutions 
can be expressed in closed form. The weakly nonlinear regime of IH convection has 
been theoretically examined in a few studies 01131II391II181 fT44l . Such analyses of 
particular solutions reveal much about bifurcations and pattern formation, but they do 
not yield robust information about heat transport. This is because the results often de¬ 
pend strongly on geometry, and also because each particular solution typically is stable 
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over only a narrow range of parameters. 

Sections 2.1 and |2.2| address each static state’s linear and nonlinear stability, re¬ 
spectively. Section |2.3| outlines a proof of f?-dependent lower bounds on the mean 
temperature for all three IH configurations. These bounds, which amount to upper 
bounds on N, are then compared with upper bounds on N in RB convection. 

The results laid out in this chapter constitute most of what can be deduced math¬ 
ematically about N and N. These results are rather meager in that they tell us neither 
the actual values assumed in the ranges 1 < < A*(R) and 1 < iV < iVfo(R), nor how 

the Prandtl number and geometry affect these values. Such questions must wait until 
the next chapter because substantial answers, so far, come only from simulations and 
laboratory experiments. 


2.1 Linear instability of static states 

In each RB and IH configuration, we can find a Rayleigh number, Rl, above which 
the static state is linearly unstable. The Prandtl number of the fluid does not affect 
this threshold. In most cases it has been proven that the linear instability is stationary, 
meaning that the non-static states that bifurcate at the point of instability are steady, 
rather than time-dependent. The method of calculating Rl is similar in every case and 
is well known from the study of the canonical RBI system. We outline this methods 
here and give references for further details. 


2.1.1 Linear stability eigenproblem 


We want to study the stability of the static state, wherein u = 0 and T = Tst{z) for the 
various Tst{z) profiles given in expression (1.12i. It is convenient to decompose the 
temperature held into its static part and a fluctuation, 0, 


7’(x,f) = 7;,(z) + 0(x,f). 

Since u and T evolve according to the Boussinesq equations ( |1.6[ )-( [TT8| l, fluctuations 
evolve according to 


V-u = 0 (2.1) 

5rU + u-Vu = —Vp-fPrV^u-|-PrR0z (2.2) 

(9f0+u-V0 = v20-r>, (2.3) 

where the prime denotes The static state enters the fluctuation equations only 
through its gradient, T/j, reflecting the fact that Boussinesq dynamics are affected only 
by relative temperature differences, not absolute temperatures. This gradient is constant 
in RB convection but varies linearly in IH convection when the heating is uniform: 

{ -1 RB1,RB2, RB3 

-z+i IHl 
-z IH2, IH3, 


(2.4) 
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where we recall that 0 < z < 1. The boundary conditions on 9 are the homogenous 
analogs of the conditions on T: 


RBI, IHl; 

0|z=O, 

0|z=l 

= 0 

(2.5) 

RB2, IH2; 

d^Q\z=l 

= 0 

(2.6) 

RB3, IH3; 

o 

II 

0|z=l 

= 0. 

(2.7) 


The fluctuation dynamics of RBI and IHl are distinguished only by differing T^f{z). 
The same is true of RB2 and IH2, and of RB3 and IH3. 


We study the stability of the zero solution of the fluctuation equations (2.1 i-( |2.3| l, 
which is equivalent to the stability of the static state. The nonlinear terms in the fluctua¬ 
tion equations can be neglected when finding linear stability thresholds. As is standard 
mi, we And a closed pair of equations governing the linear evolution of w and 9 by 


taking z • V x V x (2.2). Omitting time derivatives gives equations for the marginally 


stable states that are stationary, meaning they do not vary in time, are governed by 


VV = -RV^0 


= T>, 


( 2 . 8 ) 

(2.9) 


where := is the horizontal Laplacian operator. The validity of considering 

only stationary instabilities is discussed at the end of this subsection. 

The Rayleigh number, Rl, at which the static state becomes linearly unstable is the 
smallest R for which there is a marginally stable state—that is, the smallest R for which 


equations (2.8 l- p~9| l have a nonzero solution. This is a (generalized) eigenproblem 
whose spectrum of eigenvalues is continuous and bounded below. Assuming there 
are no horizontal boundaries, we can Fourier transform the eigenproblem in x and 
y, decomposing it into an independent eigenproblem for each horizontal wavevector 
(kxjky), where kx and ky are real. If the horizontal periods of a mode are Lx and Ly, 
then kx = 2nlLx and ky = InjLy. The resulting decomposed eigenproblems take the 
form lfT09l[T^ 




w^ > -2k^w' + k‘^w=Rk^9 
9''-k^9 = TLw, 


( 2 . 10 ) 

( 2 . 11 ) 


where v!>(z) and 9{z) are complex in general, and k^ :=k^+ky. We call k the horizontal 
wavenumber. 


The six th-order li near system (2.10l-(2.11 1 requires six boundary conditions. The 
conditions (2.51-(2.7 1 on 9 apply also to 9. The first two vi> conditions are that vi> = 0 
at both boundaries, and the other two depend on whether each boundary is no-slip or 
free-slip: 


no-slip; v!>'(0), v!>'(l) =0 (2.12) 

free-slip top; m>^(0), v!>"(1)=0 (2.13) 

free-slip bottom; vi>”(0), vi>^(l) =0 (2.14) 

free-slip; w"(0), w"(l) = 0. (2.15) 
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Throughout this work, we give results for all four pairs of velocity conditions when 
possible. In some cases, analytical bounds and experimental results are available only 
for no-slip boundaries, which are the most natural in the laboratory. Condition ( |2.13| l 
is experimentally realizable in a container with an open top. Condition ( |2.14| l might 
appear unrealizable since it describes a container with an open bottom, but it also de¬ 
scribes dynamically equivalent systems with an open top. For instance, IH convection 
with an open bottom, when viewed upside down, has the same dynamics as internally 
cooled convection with an o pen to p. _ 

For each equations (2.10i-(2.11 1 and their boundary conditions form a linear 
eigenproblem in R with a discrete spectrum that is easier to compute than the con¬ 
tinuous spectrum of (2.8 l-(|2.9|l. Th e 7?^ at which the static state loses stability is the 
smallest eigenvalue of (2.101-(2.11 1 , minimized over all admissible k^. If all horizontal 


wavenumbers are possible. 


Rl= inf R^^Uk), 
k ^>0 


(2.16) 


where 

R^^^ik) := min {R I ( |2.10| i-( |nT| ) has a nonzero solution} . (2.17) 

The definition of R^ requires an infimum rather than a minimum because the infimum 
sometime occurs in the limit k^ 0, in which case no minimum is achieved. Pertur¬ 
bations with k = 0 are not admissible since a horizontally uniform vi> would violate 
incompressibility. The value (or limit) of k at which R]^ occurs is called the critical 
wavenumber of linear instability, k^. 

Because the eigenproblem p.8| l-( [2^ is derived assuming a stationary instability, 
the resulting Rl is the value at which a steady state bifurcates from the static one. In 
cases where it is proven that all marginally stable states are indeed stationary, R > R^ 
is not only sufficient but necessary for instability of the static state. Stationarity has 
been proven for all RB cases 110211 . In the IFI3 case it follows for free-slip boundaries 
from an argument of Spiegel (see footnote 4 of M146I 1 and for no-slip boundaries from 
a theorem of Herron ll58l . The latter method of proof may suffice to show stationarity 
in the remaining IH configurations. Until that is done, we can say in those cases only 
that R> Ri is sufficient for instability. 


2.1.2 Solutions of the linear stability eigenproblem 


In RB convection, where T'^ = —1, the eigenfunctions solving (2.10l-(2.11 1 are com¬ 
binations of trigonometric and hyperbolic functions. The minimum eigenvalue at a 
given wavenumber, must satisfy an expression involving the eigenfunctions. 

This expression can be solved for R^^^{k), analytically in a few cases and numerically 
in the others 111091 fTSll . In IH convection, where varies linearly with z, the anal¬ 
ogous approach would involve hypergeometric functions, so it is simpler to solve the 
eigenproblem ( |2.10| i-( [2~TT] ) numerically. We have done this for all six configurations 
by the general method described in il4ll : discretizing the differential operators using 
a spectral collocation method and computing the spectra of the resulting matrices. Our 
computed values of Ri agree with or add precision to the values in the literature. As 
explained shortly, the exact values of Rl in the RB2 and IH2 cases can be calculated 
also by asymptotic expansion 
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Figure 2.2: Schematic diagrams of how the first marginally stable eigenvalue, 
depends on the horizontal wavenumber, k, in (a) the RBI, RB3, IHl, and IH3 cases, 
and in (b) the RB2 and IH2 cases, where heat fluxes are fixed at both boundaries. 


The {k) curve can assume one of two qual itativ e shapes in our models, depend¬ 
ing on the thermal boundary conditions. Figure 2.2 1 shows what the R^^^{k) curves 
look like in the four cases where the temperature is fixed at one or both boundaries: 
rW 

is min imized at a finite k and grows unboundedly both as k ^ 0 and as k —> oo. 
Figure [ 2 !^ shows what the R***) (k) curves look like when the temperature flux is fixed 
at both boundaries: R^^'^ approaches its inflmum as k —^ 0 and grows unboundedly as 
k —^ 00 . 

For all of our RB and IH configurations and all four pairs of velocity conditions 
( 2.12 | i-( 2T5) |, table 2.2 gives the smallest Rayleigh number, Rl, at which the static state 
undergoes a stationary, horizontally periodic instability, along with the instability’s 
horizontal wavenumber, kj^. Most of these values have been known for a long time. The 
RBI case was analyzed first, with free-slip boundaries in Rayleigh’s seminal analysis 
of 1916 Ml091 and then with other velocity conditions 11711 l^fl 021 . Later, values of 
Ri for some of our velocity conditions were reported for RB2 and RB3 II125I . IHl 
lIllllMlinaillSlIia, IH3 miiniissilsa, and then IH2 Il67ll4^. 

In most configurations, Rl is smallest when both boundaries are free-slip, largest 
when both boundaries are no-slip, and somewhere in between when one boundary is 
free-slip and the other is no-slip. The IHl configuration provides a surprising excep¬ 
tion: Ri is smallest when only the top is free-slip and largest when only the bottom is 
free-slip. 

In principle, knowing Rl is relevant to heat transport because R > R^ suggests 
that sustained convection will occur. In confined geometries, however, not all k are 
admitted, and the flow can be only approximately periodic in the horizontal directions. 
These effects raise Rl by an amount particular to the confining geometry, and only 
when R exceeds this larger value is convection guaranteed. 

The RB2 and IH2 configurations are special in that analytical expressions for Rl 
can be found for any velocity boundary conditions. This is because the critical wavenum¬ 
ber is zero (cf. table|2.2|l, so Rl can be calculated by long-wavelength asymptotics. 
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Rl 

ki 

RBI 

no-slip 

1707.76 

3.1163 

free-slip top 

1100.65 

2.6823 

free-slip bottom 

1100.65 

2.6823 

free-slip 

657.511 

2.2214 

RB2 

no-slip 

720 

0 

free-slip top 

320 

0 

free-slip bottom 

320 

0 

free-slip 

120 

0 

RB3 

no-slip 

1295.78 

2.5519 

free-slip top 

816.744 

2.2147 

free-slip bottom 

668.998 

2.0856 

free-slip 

384.693 

1.7576 

IHl 

no-slip 

37 325.2 

3.9989 

free-slip top 

16669.8 

3.0131 

free-slip bottom 

37949.4 

4.0867 

free-slip 

16992.2 

3.0277 

IH2 

no-slip 

1440 

0 

free-slip top 

576 

0 

free-slip bottom 

720 

0 

free-slip 

240 

0 

IH3 

no-slip 

2772.27 

2.6293 

free-slip top 

1612.62 

2.2611 

free-slip bottom 

1650.55 

2.1429 

free-slip 

867.766 

1.7897 


Table 2.2: Rayleigh number, Rl, above which each configuration’s static state is lin¬ 
early unstable, and the horizontal wavenumber, k^, of the linear perturbation that is 
marginally stable at R^. The RB2 and IH2 values are exact (cf. §2.1.3 1 , while the other 
values are numerical approximations that are accurate to the precision shown. 
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2.1.3 Long-wavelength asymptotics for RB2 and IH2 

In both RB2 and IH2, the infimum in the definition ( |2.16| ) of Rl occurs when A: —> 0, so 

RL = \imR^°\k). (2.18) 

k^O 

It has apparently not been proven analytically that fixed-flux boundary conditions im¬ 
ply ki — 0, so the finding must be verified on a case-by-case basis by numerically 
computing R^^\k). This has been done for RB2 and I H2, so Ri can be found exactly 
by expanding the linear stability eigenproblem (2.10i-(2.11 1 in small k^. (We cannot 
simply set A: = 0 because the limit is singular.) The eigenmode with eigenvalue (A:) 
has scaling (9(0) = k^O{w) when A: <C 1 1201 . We thus let vi> = k^W and seek solutions 
where W and 0 are both (9(1). The rescaled eigenproblem is 


= Re+2k^W" -k^W 
e” = k^ {§ + T^,w). 


(2.19) 

( 2 . 20 ) 


The no-flux conditions on 0 require the vertical integral of 0" to vanish, and this fur¬ 


nishes a consistency condition on the righthand side of (2.201: 


f\e + T:,w)dz = Q. 
Jo 


( 2 . 21 ) 


Equations ( 2.19[ )-( [Z2T| ) suffice to determine Ri, but it is possible to also retain the 
nonlinear terms in a long-wavelength expansion of the fluctuation equations. This has 
been carried out for RB2 ll20l . IH2 ll67ll . and some other convective models with fixed 
boundary fluxes Eiiia. The simpler linear calculation we describe here is contained 
in these nonlin ear analyses . 

Equations (2.191-(2.20 1 are solved asymptotically by expanding in k^: 


W{z)=Woiz)+k^W2{z) + --- 

(2.22) 

e(z) = 9o(z) +k^e2(z) + ■ ■ ■ 

(2.23) 

R = Ri k^R2 -(-■■■. 

(2.24) 


The R expansion anticipates that Rq = R^ —in other words, that A:^ = 0. We and others 
have confirmed this for RB2 and IH2 by computing the R^^^ ( k ) curv es numerically. 

The procedure for asymptotically solving equations (2.19 1 -( 2.201 to any order A:^” is 
entirely systematic. Assuming all lowe r-order terms are known, the polynomial 02„ (z) 
is found by integrating equati on (2.20 1 at 0{k^"), then the polynomial W 2 n{z) is found 
by integra ting e quation (2.19 1 at 0(P"), and finally R 2 n is found from the consistency 
condition (2.211 at 0{k^’'). Since Ri = Rq here, we only need to carry out these three 
steps at leading order to find R^. This has been done for RB2 in EOll and for IH2 in 
|. The results of the three steps are that 


Rl = 


-1 


/j T:,p{z)dz 



RB2 


IH2, 


(2.25) 


























CHAPTER 2. STABILITIES AND BOUNDS 


33 


where P{z) is the unique fourth-order polynomial that has a leading coefficient of 1/24 
and satisfies the v!> boundary conditions. For our domain of 0 < z < 1, these P{z) are 
given in ll48l . for instance. The values of Rl for various velocity conditions appear in 
table l2.2l above. 


2.2 Energy stability of static states 

Knowing the Rayleigh number, R^, above which a static state is linearly unstable would 
be well complemented by knowing the critical Rayleigh number, Rc, below which it 
is the globally attracting state of the fully nonlinear dynamics. This is difficult in 
general, so we settle for finding a so-called energy Rayleigh number. Re, that is a lower 
bound on Rc- Since linear instability implies nonlinear instability, we can anticipate that 
Re Rc ^ Rl- RB convection is special in that Re = Rc = Rl 1741 . IH convection is 
more complicated in that Re < Rl for the largest known Re, as depicted in figure [ZT] 
The values of Re and Rl serve as upper and lower bounds on Rc, respectively, that hold 
uniformly for all Pr. The exact values of Rc can depend on Pr and are not yet known. 
The upper bounds on Rc can be tightened by finding particular subcritical solutions, 
as has been done for IHl EH and IH3 lfTT8lfT44l . since any R at which subcritical 
convection persists must be larger than Rc- Proving a lower bound tighter than Re is a 
more daunting challenge. 


2.2.1 Lyapunov stability and the energy method 

The global stability of the static state is equivalent to the global stability of the zero 
solution to the fluctuation equations (|ZT]l-(|^. The nonlinear terms in those equations 
that could be ignored in the linear stability analysis must now be included. The typical 
method of proving global stability, due to Lyapunov, requires finding a functional of 
the state variables that is positive definite and whose evolution is negative definite. That 
is, we must find a functional .if [u, 0] such that 

.if[u,0]>O (2.26) 

^.if[u,0]<O. (2.27) 


To show also that the static state attracts all initial conditions, it suffices for the above 
inequalities to be strict whenever u or 0 is nonzero. In the convective systems we are 
studying, the best we can hope for is to find an .if where (2.26 1 -(2.27 1 hold for R below 
some finite value, R That is. 


R^ '■= sup{R I .if satisfies ( |2.26| l-( |ZZ7l )} . (2.28) 


The critical Rayleigh number, Rc, is the largest R at which any Lyapunov functional 
exists, 

Rc := sup R . (2.29) 


There is no universally successful method for constructing Lyapunov functionals, 
let alone the optimal .if that is valid for R up to Rc- It is even difficult to confirm that 
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an optimal ^ is indeed optimal, except when Rc = Rl, as in RB convection. All we 
can do in general is make educated guesses for determine the corresponding values 
of and declare the largest R^ we can find to be a lower bound on Rc- In fluid 
dynamical systems like ours, even this guess-and-check procedure cannot be carried 
out for general .jSf because it is too difficult to determine whether the second Lyapunov 
condition ( |2.27| i holds. In most nonlinear analyses of fluid stability, this trouble is 
avoided by considering only a particular subset of possible Lyapunov functionals for 
which it is tractable to check the second Lyapunov condition. This approach is called 
the energy method. 

The energy method in fluid mechanics II 1201 17^ 11331 is a special case of Lya¬ 
punov’s method. In one definition of the energy method that is neither the narrowest 
nor the broadest definition possible, the Lyapunov functional, which is called the en¬ 
ergy, has two special features: 


1. The energy is quadratic in the state variables. 


2. The energy is conserved by the nonlinear terms of the fluctuation equations (2.2 1 - 
(2.31, meaning that these terms do not contribute to the expression for the time- 
evolution of the energy. 


The energy method is so named because quadratic quantities are often proportional to 
physical energies. Here we follow Joseph fMl in considering energies of the form 


Ey[u,9]{t) i 


(prR 


u|2-P702) 


dx. 


(2.30) 


where / denotes an instantaneous volume average. The constant 7 > 0 is called a 
coupling parameter, and each positive value defines an energy that is a valid Lyapunov 
functional for R up to some This value of is maximized by some optimal 
choice of 7 , where it achieves the critical Rayleigh number of energy stability. Re'. 

Re '.= maxsup {R | Ey satisfies ( |2.26| l-( |2!27] i| . (2.31) 

The value of Re is the best lower bound on Rc that we find by the energy method, 
though it is likely still smaller than Rc. Deriving a better lower bound on Rc would 
require going beyond the energy method to search over a larger class of Lyapunov 
functionals, and this presents technical challenges. Progress beyond the energy method 
has been made for a few shear flow models 117811231 but not yet for a convective system. 


2.2.2 Energy stability eigenproblem 


The functional Ey suffices to show that the static state is globally stable whenever it 
satisfies conditions ( 2.26[ )-( [2!27] l. The first condition holds whenever all the parameters 
are positive, so it remains only to determine the parameters for which iEy< 0. Adding 
the volume averages of u • ( 2.21 and yd x (2.3) and then integrating by parts gives 


TtEy = -j [^1 Vu|2 + 7|V0|2 - (1 - 77;;) we] dx. 


( 2 . 32 ) 
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The static state is globally attracting if the righthand side of { 232) is negative definite. 
Like the linear stability threshold, the satisfaction of this condition depends on R but not 
on Pr. Only static states have this feature; other solutions and their stabilities depend 
also on Pr. 

The calculus of variations yields a necessary and sufficient condition for the right- 
hand side of (2.32 1 to be negative definite. In particular, Ey is a Lyapunov func¬ 
tional when R is smaller than all eigenvalues, R, of the (generalized) eigenproblem 

ImlEllml 


- 2k^w” + k‘^w = \Rk^ (1 - yr,;) 0 (2.33) 

7(0"-k20) =-i(l-77;;)M). (2.34) 

The boundary conditions are the same as in the linear stability eigenproblem of §2. L2[ 
and again vi>(z) and Q{z) can be complex, and k is the horizontal wavenumber. Expres¬ 
sion ( |2.31| l for Re can thus be restated as 

Re = max inf min {R | ( |2.33| l-( [234l l has a nonzero solution} . (2.35) 

It is a special feature of the energy method, and not of Lyapunov’s method in general, 
that the nonlinear stability analysis can be reduced to the solution of a linear eigen¬ 
problem, much like the linear stability analysis. 


2.2.3 Solutions of the energy stability eigenproblem 

In RB convection, where T/, = — 1, there is no need to solve the energy stability eigen¬ 
problem ( |2.33| l-(2.34 1 because it is identical to the linear stability eigenproblem ( 2.10| i- 
( |2.1 l| i, so long as the energy is defined with 7=1. This energy is thus a valid Lya¬ 
punov functional for all R up to Re- (The agreement of the two eigenproblems reflects 
the symmetry of the linear stability operator; see ll45lI133L 1 We expect 7=1 to be 
the optimal coupling parameter since Re should not exceed Re, and indeed this can be 
shown directly EH. These observations justify our earlier assertion that Re = Rc= Rl 
in RB convection, making subcritical instability impossible. 

In IH convection. Re must be calculated by performing the double optimization 
of expression ( 2.35| l, which requires solving the eigenproblem (2.331-(2.34 1 . In all 
IH cases the strict inequality Re < Re holds. Table [23] gives values of Re for the 
various IH configurations, along with the percent differences between Re and Re, and 


the arguments, ']4 and kE, that yield the maxima and infima in expression (2.35 1 . 

The relative magnitudes of the gaps between Re and Re are on the order of 1 % in 
IH2 and IH3, where the bottom is insulating, but are much larger in IHl, where heat 
escapes across both boundaries. We cannot say whether the larger gaps in IHl are 
necessitated by subcritical solutions or are only mathematical artifacts of the optimal 
energies being poor approximations of the truly optimal Lyapunov functionals. The 
most energy-unstable wavenumber, kg, is fairly close to in IHl and IH3, and k^ = 
k^ = 0 in IH2. The optimal coupling parameters, y*, are all significantly larger than 
unity, which is their optimal value in the RB cases. 
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Re 

% below Re 

f 

kn 

IHl 

no-slip 

26 926.6 

27.9 % 

8.8831 

3.6174 

free-slip top 

12 620.2 

24.3 % 

7.9626 

2.9014 

free-slip bottom 

24722.8 

34.9 % 

9.1975 

3.3664 

free-slip 

10618.1 

37.5 % 

8.8516 

2.5498 

IH2 

no-slip 

1429.86 

0.704 % 

1.9720 

0 

free-slip top 

573.391 

0.453 % 

1.7838 

0 

free-slip bottom 

714.929 

0.704 % 

2.2185 

0 

free-slip 

239.055 

0.394% 

1.9843 

0 

IH3 

no-slip 

2737.16 

1.27% 

2.0678 

2.6355 

free-slip top 

1594.42 

1.13% 

1.9185 

2.2661 

free-slip bottom 

1624.26 

1.59% 

2.3702 

2.1512 

free-slip 

855.674 

1.39% 

2.1821 

1.7958 


Table 2.3: Rayleigh number (Re) below which the energy method proves that each 
IH configuration’s static state is globally attracting, the percentage of Re by which Re 
falls short of Re, the optimal coupling parameter ( 7 *) used to define the energy that 
is a valid Lyapunov functional for all R < Re, and the horizontal wavenumber (A:^) 
at which the infimum in (2.35i occurs for the optimal energy. The IH2 values are 
numerical approximations to the exact analytical expressions ( |2.36| l. The IHl and IH3 
values are computed numerically and are accurate to the precision shown. 
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For the IHl and IH3 cases, we numerically computed the values of Re, Y*, and 
kE that appear in table 2.3 using the same spectral collocation method that we used to 
compute Re- A similar energy stability analysis was carried out by Kulacki & Goldstein 
1861 . The values of Re that they reported are smaller than our own, perhaps because 
their coupling parameters were not quite optimalj^ An energy stability analysis was 
carried out more recently for the IH3 configuration with no-slip boundaries M132I . and 
those findings agree exactly with our own. 

In the IH2 case, the energy stability eigenproblem can be solved exactly using 
long-wavelength asymptotics, much like the linear sta bility eigenproblem (cf. §2.1.3 i. 
This is possible because the infimum of expression (2.35 1 is reached as 0, an 
observation that has not been proven but has been confirmed numerically ll4^ . The 
asymptotic calculations, which are detailed in HSl . give the exact expressions 


Re 


'2880(6^35-35) 

360 ( 9 ^^- 175 ) 

1440(6^35-35) 

1440(8^7-21) 


no-slip 
free-slip top 
free-slip bottom 
free-slip 


(2.36) 


Numerical approximations of the above values appear in table |2. 3 [ above. 

From the standpoint of scientific and engineering applications, the value of know¬ 
ing Re in IH convection is that we know convection cannot be sustained when R < Re- 
When R lies between Re and Re, little is known about when convection can occur, 
apart from some instances of subcritical convection that have been computed in IHl 
11431 and IH3 0118111441 . This ambiguous regime between Re and Re is small in IH2 
and IH3, and thus of not much practical importance, but it is much larger in IHl. In 
any event, we cannot claim to fully understand the static states until we know when 
subcritical convection is possible—that is, until we know the true value of Rc for every 
Pr. Lower bounds on Rc could be improved by looking beyond the energy method to 
find better Lyapunov functionals, and upper bounds could be improved by numerically 
computing steady states that exist in the subcritical regimes. 


2.3 Bounds depending on the Rayleigh number 

Our main goal is to predict the parameter-dependence of integral quantities like (wT), 



they give necessary and sufficient conditions for N and N to equal unity. In particular, 
R < Re guarantees that both quantities equal unity, and R > Re guarantees that both 
are greater than unity. At large R, where convection is strong and complicated, exact 


*What we call the IHl case is designated in (86) by the parameter Biq — co, and what we call the IH3 
case is designated by the parameters Biq — 0 and Bii — 0 °. Their Rayleigh numbers are converted to our 
scaling upon multiplication by 64. 
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expressions for N{R,Pr) and N{R,Pr) are not available. Instead, we seek to bound 
these quantities analytically. 

The only parameter-dependent bounds that have been proven for RB or IH config¬ 
urations can be stated as upper bounds on how quickly N or N, respectively, can grow 
as R is raised. Upper bounds on N have not been proven in IH convection but seem 
likely to hold (cf. j |1.6.4| i. We cannot improve the lower bounds of unity since known 
techniques cannot distinguish realizable solutions from the unstable static states. 

In this section we outline a proof of lower bounds on the mean temperature, 5(7’) := 
(r — Tt), in IH convection. (Recall that 5(7’) = (T) in IHl and IH3 but not in IH2, 
and that lower bounds on 5(7’) are equivalent to upper bounds on N.) Our exposition 
combines existing results for IHl 1921 and IH2 l4^ and a new result for IH3. The proof 
employs the background method l3^l32l . which requires no assumptions beyond the 
governing equations. Like similar variational methods l62l flOl l63l . the background 
method makes progress by relaxing the constraints on u and T. Instead of enforcing 
the full Boussinesq equations, we enforce only incompressibility, the boundary con¬ 
ditions, and a few integral relations that follow from the governing equations. This 
yields bounds that hold for an enlarged class of u and T that includes solutions of the 
Boussinesq equations. 

Two main integral relations are typically enforced when the background method is 
applied to convection. They are called the power integrals and are derived by taking 
(u • (1.71) and {T x (1.8)) and integrate by parts to find 


(|Vu|2) =7;(w7’) 

{ l-l-(wr) RBI 

5r = 1 - (wT) RB2, RB3 
5(7’) IH1,IH2, IH3. 


(2.37) 

(2.38) 


Time derivatives have vanished from the above relations in the infinite-time limit since 
the volume integrals of |u| and |7’| are bounded uniformly in time. This boundedness is 
proven as a byproduct of the background-method analysis itself ll37ll^ . The absence 
of Pr from the relaxed constraints on u and T precludes our analysis from producing 
bounds that depend on Pr. 

In all six RB and IH configurations, the bounds that have been proven by the back¬ 
ground method amount to bounds on the thermal dissipation, ^ | VTp), though they are 
often stated in terms of quantities like {wT) or 5(7’) that are related to (|V7’|^) by 
( 2.38| l. The thermal dissipation is bounded above in RBI and below in the other five 
cases. These results constitute upper bounds on N in RB convection and upper bounds 
on N in IH convection. 


2.3.1 Proof by the background method 

We now prove for all three IH configurations that the dimensionless mean temperature, 
5(7’), decays no faster than R^^/^. This is equivalent to the dimensional mean temper¬ 
ature, 5(7’) A, growing with the rate of volumetric heating, 77, no slower than 77^/^. We 
assume a no-slip top in the IH2 case but need not do so in the IHl or IH3 cases. 
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To apply the background method, we decompose the temperature into a so-called 
background profile, t(z), and the remaining part, 0(x,f): 

T{x,t) = T:{z) + &{x,t). (2.39) 


The bound we obtain depends on the t(z) we choose. The background profile does not 
generally solve the governing equations, in which case 0 does not evolve according to 
the fluctuation equations ( |2.1| i-( [Z^ . 

The t(z) we choose must satisfy several conditions. First, it must be continu¬ 
ous. Second, it must satisfy the same boundary conditions as T since this lets 0 sat¬ 
isfy the corresponding homogenous conditions. In practice, however, only the fixed- 
temperature conditions on t(z) need to be enforced. This is because fixed-flux condi¬ 
tions on t(z) can be met by boundary layers whose influence vanishes as we send their 
thicknesses to zero ll48l . These limiting bounds are the same as those reached by sim¬ 
ply ignoring the fixed-flux conditions on t(z), so we do the latter in our calculations. 
Finally, t{z) must be chosen to make a particular quantity nonnegative, as explained 
below. We will see that for all admissible t{z), 

5{T)>2{t-Xt)-{x'^) ■ (2.40) 


We choose simple t(z) that make our calculations analytically tractable, thereby yield¬ 
ing analytical bounds that are valid at all R. Optimizing t(z) numerically at a given R 
would give a tighter bound (as in M103I 1. but the bound would apply only at that value 
oIR. 

To see where the inequality (2.40|l comes from, and when it holds, we expand the 
power integral (|2.38|l for the IH cases using the decomposition (|2.39|l to find 


8{T) = (ivrl^) = (t'2) +2(t'0') + (|V0|2) , 


(2.41) 


where primes denote z-derivatives. Our goal is to bound 5(7’) below. (We could equally 
well speak of bounding (| VTp) below or bounding N above.) The (t'0') term in the 
above expression is difficult to bound, so we eliminate it using a third and final inte¬ 


gral relation. Integrating t(z) against the temperature equation (1.8i gives the needed 
relation ll92l . 

(t'0') = (t - Tt) - (t'^) -h (t'w0) , (2.42) 

where the top temperature Xt may be nonzero only in the IFI2 case. Eliminating (t'0') 
from expression ( |2.41[ ) gives 

5(7) =2(T-Tr)-(T'2)-h(|V0|2)-h2(T'w0). (2.43) 


From the above equality it follows that the lower bound p.40| l would hold if we could 
show (|V0p) + 2{ t'w0) > 0. This is an impossible task for arbitrary w, however, 
since the velocity enters only in the sign-indefinite term. Apparently, the temperature 


power integral (2.381 alone is not suf ficient ly constraining. We need the additional con¬ 
straint of the velocity power integral (2.37 1 , which tells us that a (| Vu p) — (wT)) = 
0 for any a. Adding this relation to ( 2.43[ ) shows that the lower bound ( 2.40| i on 5(7’) 
would follow from the nonnegativity of the quadratic functional 

^[u,0;t(z,7;)] := |(|Vu|2) + (|V0|2) + ((2t' -a)w0) . 


(2.44) 
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Figure 2.3: Schematic of the class of background profiles, t(z), that we consider. The 
parameters a, b, and 5 are optimized, within some constraints, to maximize the lower 
bounds on the mean temperature. 


We must choose a t(z) for which we can verify that ^ > 0 for all admissible u and 0. 

More generally, the background method is carried out by finding an expression, 
equal to the quantity to be bounded, that takes the form where ^ is a functional 

of the the background field alone, while ^ depends also on the other fields. The key 
idea is that ^ will be a lower bound when we can show that ^ is nonnegative (or 
an upper bound when we can show that ^ is nonpo sitive). In the present analysis, 
^ := 2 {t — Tt) — and ^ is as defined in (2.44i. 


Two objectives compete in the choice of t(z): making the lower bound (2.40i as 
large as possible, and maintaining the nonnegativity of ^ that is needed for that bound 
to be valid. Here, we optimize t(z) only among profiles consisting of two linear pieces. 
Such profiles can all be written in the ansatz 


0<z<l-5, 




(2.45) 


where the geometric meanings of parameters a, b, and 5 are shown in figure |2.3[ We 
will see that the top piece of t(z) is a boundary layer whose thickness, 5, goes to 
zero as R ^ oo. The bottom piece of t(z) has a slope that is half the value of the yet- 
unspecified constant a, a known trick Il32ll92l for making the sign-indefinite term of ^ 
vanish outside the boundary layer. 

For our three-parameter family of background profiles p.45 1, the lower bound 
( |2.40| i becomes 

5{T) >b{2-5) + ‘i{l-5)-[i+ab) (i - 1) - f. (2.46) 

With a no-slip top in IH2 and any velocity conditions in IHl or IH3, it can be shown 
that ^ > 0 is satisfied when 5 is no larger than | 


d^ = 


\64aR-^ IH1,IH3 
1 32aR-^ IH2. 


(2.47) 
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We choose this 5 because the tightest bounds result from choosing 5 as large as possi¬ 
ble. 

In the IH2 and IH3 cases, we are free to choose the parameters 5, a, and b to max¬ 
imize the lower bound ( |2.46| l subject to \2A1) . In IHl, the lower boundary condition 
requires that b = 0, so we are free to choose only 5 and a. This maximization is car¬ 
ried out for IHl and IH2 in ll92l and HSl . respectively, and the procedure for IH3 is 
analogous. The resulting optimal parameters are 


1 

f 4R-'/3 

r 4R-V3 i 

5* = 

12'/3r-i/3 a* = 1 

3.12>/3^_1/3 _ I 

1 


i ^7;->/3 1 


0 


^^- 1/3 


^^- 1/3 


IHl 

IH2 

IH3, 

(2.48) 


for which the lower bound (|2.46|l becomes 


d{T)> 


R-U^- AR-^!^ IHl 

(3)V=„-./3_§(3)W„-V, 

3^^-l/3_ 89|(i^-2/3 ih3. 


(2.49) 


At large R, the leading-order terms of the bounds dominate; 


5(7') > 


R-V3 

IHl 


l.2SR-^/^ 

IH2 

(2.50) 

1.62R-^/3 

IH3. 



When Lu et al. Il92l proved the above bound for IHl, they also raised the prefactor 
from 1 to 1.09 by generalizing the ansatz of t(z) to include a bottom boundary layer, 
although this required solving an algebraic equation numerically. Their proof carries 
through for the IH3 case also, so they in fact proved the asymptotic lower bound of 
1.09 7?^ for both IHl and IH3. By dropping the condition t( 0) = 0 in IH3, where it 
is not needed, we have raised the prefactor to 1.62. Optimizing t(z) beyond our limited 
ansatz would lower the prefactors of the bounds, but results of numerically optimizing 
t(z) in the RBI case suggest that the scaling of the bounds would not change 01031 . 


2.3.2 Similarities between RB and IH bounds 


A main virtue of the way we have defined the Nusselt numbers N and N and the diag¬ 
nostic Rayleigh number Ra and Ra is that bounds for the various configurations all have 
the sa me scaling when expressed using these quantities. Recalling that the definitions 
(1.33 1 of TV are inversely proportional to 5(7’), and that Ra = R/N in IH convection. 


we see that the asymptotic bounds (2.50 1 on 5(7) become 


N< 


0.025 

IHl 

0.132^'^^ 

IH2 

0.094 Ra ' 

IH3 


(2.51) 
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These upper bounds on N appear in table 2.1 at the start of this chapter, along with 


the best known bounds on N in the RB configurations. The RB bounds have different 
prefactors but the same exponent, scaling proportionally to The RBI prefactor 

in table |2.1 [ comes from the improvement on ll^ by Blasting & Kerswell 11031 . who 
also showed that their bound could not be improved without additional constraints. The 
prefactor in the other two RB cases comes from 1991 . where the analysis was aimed at 
RB2 but carries through for RB3 also. 

Upper bounds with an exponent of 1 /2 are the best available for three-dimensional 
convection in general, but bounds with smaller exponents have been proven in special 
cases. Here too, analogies hold between various configurations if results are stated 
in terms of our diagnostic parameters. When the boundaries are free-slip, and either 
Pr = oo or the flow is two-dimensional, upper bounds with exponents of 5/12 have been 
proven for RBI fll50lll5l1l and IHl 1149111511 . When Pr = oo with no-slip boundaries, 
the best known bounds on N scale like Ra*/^(loglogRa)'/^ in RBI 1 1001 and like 
Ra*/3(logRa)*/2 in RB2 and RB3 IB^ . and the best known bound on N scale like 

Ra ' (logRfl)*/^ in IHl 11491 . Bounds with exponents smaller than 1 /2 are yet to be 
reported for the other RB or IH configurations. 

Now that we have seen how the background method works, we can understand 
why it is challenging in the IH cases to prove upper boun ds on {wT). This quantity 
is related to (|Vup) by the velocity power integral (2.37 1 but is not is not related a 
priori to (iVTp) in IH convection. However, the parameter-dependent bounds that 
have been proven for convective models all amount to bounds on (|Vrp) and rely 
on a background decomposition of the temperature held. Bounding /|Vup) instead 
suggests a background decomposition of the velocity held, which has been carried out 
for shear flows (e.g. in Ill]) but not for convection. 
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The first two chapters have summarized features of heat transport in IH convection 
and RB convection that can be ascertained analytically from the Boussinesq equations. 
In this final chapter we summarize findings on IH convection from physical and com¬ 
putational experiments. Analogous results for RB convection are described only min¬ 
imally, as the experimental literature on RB convection is vast and has been reviewed 
elsewhere (e.g. lfT22ll47imi^l28lli. 

Precise laboratory experiments on IH convections are inherently more difficult to 
carry out than similar experiments on RB convection. Both require maintaining the 
chosen thermal boundary conditions, but IH experiments also require producing heat 
internally in a controlled way. If our simple models are to apply, the heat production 
should be constant and uniform. In most experiments, the internal heating has been 
achieved by Joule heating, where the working fluid is an electrolytic solution that is 
heated by passing current through it. Two sets of experiments P107ll88l used a different 
method, wherein heating elements were distributed throughout the domain. Although 
neither method heats uniformly, it is possible that rapid mixing by strong convection 
limits the influence of non-uniformity. This is supported by the fairly good agreement 
between non-uniformly heated experiments and uniformly heated simulations. 

Numerical simulations of IH convection avoid unknown variations in heating rate or 
material properties. However, most numerical studies were carried out several decades 
ago and were limited to 2D and fairly small R. The larger values of R accessible with 
modern computers have been simulated only a few times, and much of the parameter 
space that could now be reached has yet to be explored. 

Experimental findings before 1985 are collected in the review of Kulacki & Richards 
187], who discuss findings on IHl, IH3, and some similar configurations. The slightly 
later review of Cheung & Chawla ESI adds various scaling arguments for heat trans¬ 
port. Nourgaliev et al. 19^ summarize heat fluxes in these same early experiments, as 
well as in experiments with curved geometries and cooled side walls. 

A number of experiments have examined heat transport quantitatively, and a num¬ 
ber of others have focused on qualitative pattern formation near the onset of convection. 
Here we cite studies of both types but focus on quantitative findings, and we restrict 
ourselves to experiments that closely resemble one of the IH configurations defined in 
figure |1.1| Convection with internal heating has been studied also with various com¬ 
plications that we do not confront, such as cooled side walls 01281 [H |5] [35] |3 [H 

I^ll21ll^l22ll44l . non-uniform heating 01 1 Illl32lll38l[82] . self-gravitating spheres 
0112111 14lfm[T3]l6^ll 101 . and hybrid configurations driven both internally and by the 
boundary conditions ifT^ iTTl l75l l95l IT^ fT^ IM1 l2l fT^ [T45l l55l l9l 1^. 


Section |3T| addresses IH3, the last of the three IH configurations defined in figure 
O Section]^] addresses IHl, which is in some ways more complicated than IH3. 
We are not aware of any heat transport findings on the IH2 configuration, though 2D 
simulations have been carried out to study scale selection EllIlT]. Section [33] suggests 
directions for future work. 
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3.1 The IH3 configuration 

The internally heated configuration we call IH3, which is bounded above by a perfect 
conductor and below by a perfect insulator, has been the subject of numerous laboratory 
experiments 014211 1 1911^l43iri0811861l83lll07l[T37ll88ll 1361 . as well as computational 
studies both in 2D |[l39l|95l[39l|67l and in 3D ||T39l[ffil|6l][n6l|65l[l7l[l6l. Many of 
these investigations have focused on pattern formation and scale selection, which we 
do not discuss here. Our interest is in quantities relevant to heat transport, including the 
mean vertical temperature profile, T{z), and the mean temperature difference between 
the boundaries, 5T. No data are available on the mean fluid temperature, 5(7’). 

Computational studies of IH3 convection that report T{z) or 5T are all several 
decades old. Most are limited to the steady states that are stable at modest R 1951 
1139111441II 161 . At larger R, unsteady 2D simulations have been carried out using a 
turbulence closure model HOl and by direct numerical simulation (DNS) 1^ . although 
some runs in the latter study appear under-resolved. As far as we know, unsteady IH3 
convection has not been simulated in 3D. The largest R that has been reached in 2D 
DNS of IH3 could be greatly exceeded in 3D DNS on modern parallel computers. 

Laboratory experiments, most of which were carried out in the 1970s, furnish 
nearly everything we know about T{z) and 57 in IH3 convection at large R l43lll081 
l86l l8^ ri07l [88l . These findings are subject to the uncertainties inherent to IH ex¬ 
periments, so there is cause to repeat them numerically. The largest R reached in past 
laboratory experiments of IH3 could now be approached by 3D DNS, albeit in a smaller 
spatial domain. 

3.1.1 Temperature profiles 

Several authors have reported mean vertical temperature profiles. It is simple to ob¬ 
tain T{z) in numerical studies by averaging steady flows horizontally 01391 l95l 11441 
or averaging unsteady flows both horizontally and temporally ll^l40ll . In laboratory 
experiments, vertical profiles have been obtained by measuring temperatures at fixed 
points and averaging only over time 01081 . If the flow is horizontally isotropic in a 
statistical sense, and time averages are sufficiently long, then the same mean profile 
would be obtained whether or not horizontal averages are also taken. This is generally 
expected to be true at large R when side walls are absent or negligible. Some tran¬ 
sient profiles have been reported also ll86l [^ . but these do not bear directly on the 
infinite-ti me av erages we seek. 

Figure [3T| shows T (z) profiles for 2D steady states computed by Thirlby 01391 for 
Pr = 6.8 and relatively small R. As R is raised and convection strengthens, the dimen¬ 
sionless temperature decreases, and the interior becomes closer to isothermal. When 
convection is sufficiently strong, the maximum value of T (z) occurs inside the layer, 
rather than at the bottom boundary. At still larger R, where convection is stronger 
and unsteady, the experimentally measured profiles of Ralph & Roberts 01081 follow 
similar trends but are closer to isothermal in the interior, lacking the pronounced tem¬ 
perature inversion found in the steady states of figure |3.1| The unsteady motions re¬ 
sponsible for homogenizing temperature outside the thermal boundary layer are evident 
in the IH3 temperature field of figure [L^; plumes emerge from the unstably stratified 
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Figure 3.1; Numerically computed mean temperature profiles, T{z), for steady 2D 
convection between no-slip boundaries. The Prandtl number is 6.8. The rightmost 
profile is that of the static state. The others, from right to left, are for = 10^ x 
(3,4,5,10,20,30,52). The hgure is adapted from hgure 2 of Thirlby 11391 


upper boundary layer and strongly mix fluid in the rest of the domain. 


3.1.2 Mean temperature differences 

The difference by which the dimensionless temperature at any point exceeds its average 
value at the top boundary, T — Tj, tends to decrease as IH convection strengthens. This 
fact u nderlies the two related but distinct measures of convective strength discussed in 
) 1.6 5(7’), which is the average of T — Tj over time and the entire volume, and 5T, 


which is its average over time and the bottom boundary. (Recall that ST is also the 
mean vertical conduction, and in IH3 it is tied to the mean vertical convection, (wT), 
by the relation 5T + {wT) — 1 /2.) Both 5T and 5(7’) are maximal in the static state 
and most likely approach zero in convective flows as 7? —^ 0 °. 

A primary question for experimentalists is how quickly 57’ and 5(7’) fall as R is 
raised, along with how this answer is affected by the velocity boundary conditions, 
Prandtl number, and geometry. If 5r and 5{T) vary approximately as powers of R 
when other parameters are held constant, data can be captured by hts of the form 


5T ^aR- 


5{T)r^bR- 


(3.1) 


In the IH3 case, the Nusselt numbers and diagnostic Rayleigh numbers dehned in 

§ iT:64lfL63] are 


A = 1/257’ 
A= l/35(r) 


Ra=R/N 

'Ra=R/N. 


(3.2) 

(3.3) 
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Pr 

R 

5rfit 

Vfrt 

Laboratory experiments 

Fiedler & Wille US 

6-7 

o 

1 

o 

1.90R^°-528 

0.177Ra‘’-295 

Ralph & Roberts ||l08l 

6-7 

2.3-10^-6.0-lO*^ 

2.62R-‘’'25 

0.110Ra‘’-55 

Kulacki & Nagle gb) 

6.2-6.6 

1.5-105-2.5-10‘^ 

3.28R-°'539 

0.0845 

Kulacki & Emara if^ 

2.7-6.9 

1.89-10^-2.17-10'2 

2.53R-°-527 

0.123Ra‘’'294 

Ralph et al. ||107| 

6-7 

10‘^-7-10*^ 


cRa^-^^ 

Lee et al. 1(8^ 

0.71-0.74 

9.9-10‘^-3.3-10" 

2.84R-0-247 

0.0996 Ra“'528 

Simulations (2D DNS) 

Mckenzie et al. If95l 
(free-slip, steady) 

oo 

1.2-10^-7.0-105 

aR-^-^^ 


Emara & Kulacki 1^ 

(free-slip top) 

6.5 

oo 

O 

1 

o 

1.07R-0'*2 

0.397 

Emara & Kulacki ll^ 

(no-slip) 

6.5 

5-105-5-10* 

2.38R-°'223 

0.134Rfl‘’'287 

oiwi m 

(no-slip, steady) 

6.5 

o 

1 

o 

oo 

3.07R-°-255 

0.0876 Rfl“-5'‘2 


Table 3.1; Summary of IH3 experiments and simulations reporting approximate power- 
law dependence of 5T on R. Internal heating was achieved by electric curTent in the 
first four experiments and by heating elements in the last two. The Prandtl number 
range 6 — 7 is an estimate for experiments that used aqueous solutions but did not 
report Pr measurements Il43lll08lfl07l . 


Various authors have considered quantities like N in past studies of 1H3, and Fiedler 
& Wille 1431 considered both N and Ra together. Restated in terms of the diagnostic 
variables, the fits of expression (|3.1[) become 


N ~ cRa^ 


N ^ dRa , 


(3.4) 


were y = a/{I — a), d — p/{I — fi), c = and d = 

The mea n tem perature difference 5T has been measured in a number of experi¬ 
ments. Table 3.1 summarizes past fits of the form 5T ~ aR^“, along with their cor¬ 
responding re-expressions as fits of the form N ^ cRaT. Ranges of Pr and R are also 
given. The stated ranges of R are those over which the data have been fit. The Prandtl 
number would ideally be held constant as R is changed, but slight variations are un¬ 
avoidable in the laboratory. Numerical studies do not suffer from this uncertainty, but 
the simulation results in table |3. 1 [ nonetheless must be regarded with care since they all 
are 2D and seem to be somewhat under-resolved at larger R. 
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The decay rates of 5T reported for the six laboratory experiments in table 


3.1 


fall 

between a = 0.227 and a — 0.25. This means that the dimensional temperature differ¬ 
ence between the boundaries, 5TA, grows with the volumetric heating at rates between 

^o.vs^ndi/o™^ 

When the 5T fits in table 3.1 are restated in the form N ~ cRa^, the exponents 
fall between 7 = 0.294 and 7 = 0.33. This range agrees very well with the analogous 
range of 7 measured in RBI experiments, where fits still take the form N ~ cRa^, but 
with ;= 1 + {wT) and Ra := R (cf. i 1.6.4 1 . The RBI exponents summarized in table 
1 of 


lie between 0.25 and 0.33, excluding the very small Pr values for which 
corresponding IH3 data are unavailable. Exponents larger than 0.33 have sometimes 
been measured in RBI experiments at very large R 1121115^ . but no IH3 experiments 
have reached such R values. The similarity between measured values of yin RBI and 
IH3 is one of the analogies brought out by our chosen definitions of N and Ra. 

No data have been reported on the volume-averaged quantity 5(7’), so we can¬ 
not say exactly what exponents would emerge from fits of the form 5(7’) ~ 57?^^ or 
~ ~S 

N ^ dRa . We can reasonably estimate the exponents, however, since the temperature 
profiles that have been reported are close to isothermal outside their boundary layers. 
This suggests that the values of 5(7’) and 5T become ever closer as 7? 00 , in which 

case a Ki p and 7 « 5 for sufficiently large R. This speculation remains to be tested 
since volume averages like 5(7’) are difficult to measure in the laboratory. They are 
easy to extract from simulations, however, and we hope that future numerical studies 
will report 5(7’). 

Whereas we have data on 57’ but not on 5(7’)—or, equivalently, on N but not on 
N —the state of affairs for analytical bounds is just the opposite. We have conjectured 
in chapterbut have not proven, than N obeys an upper bound of the form cRa^/^. 
The experimental exponents j3 in table 3.1 are all smaller than 1/2 and thus consistent 
with this conjecture. On the other hand, we have proven in chapter chap; stab that N 


can grow no faster than 0.093 7?a 
IH3 configuration. 


1/2 


, but no data on 5(7’) have been reported for the 


3.2 The IHl configuration 

The internally heated configuration we call IHl, which is bounded above and below by 
perfect conductors of equal temperature, has been studied in the laboratory ll84llMll94l 
[MllTOllMlIll, as well as numerically both in 2 D |l 68 l El [M] [OSl [Hll EH EB and 
in 3D Il5^ll541l50ll . Almost all of these studies have reported quantitatively on heat 
transport in some way. 

Numerical computations of IHl include both steady states ifToTlfmfTrl and 
DNS. Whereas DNS of the IH3 configuration has been limited to a single 2D study, 
DNS of the IHl configuration has been carried out up to fairly large R in both 2D 
IMJESHm and 3D ll5^IT541 1501 . 
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3.2.1 Temperature profiles 

Mean vertical temperature profiles have been reported in a number of studies. Numeri¬ 
cal studies provide profiles, T{z), that are averaged horizontally and, if the simulations 
are unsteady, over time as well OlOll l94l 11351 [39l l53l 11541 ISTl l50ll . In the labora¬ 
tory, profiles measured pointwise by temperature probes are averaged only over time 
O107l[88l . while profiles gleaned from interferograms are instantaneous but effectively 
averaged over a horizontal direction llMll94l . 

Figure [ 3 !^ shows T (z) profiles from the 3D DNS data of Goluskin & van der Poel 
iBOl . As in the IH3 configuration, raising R strengthens convection, which decreases 
the dimensionless temperature and brings the interior closer to isothermality. When R 
is large enough for thermal boundary layers to be discernible, the top boundary layer 
is visibly thinner than the bottom one. This reflects the up-down asymmetry of heat 
fluxes; more of the produced heat flows outward across the top boundary than across 
the bottom one, as quantified in the next subsection. The same basic features are evi¬ 
dent in Figure [T2] 3, which shows an interferogram from the experiments of Kulacki & 
Goldstein IIMl . The interferogram measures horizontally averaged optical properties of 
the fluid that vary with its temperature, and a line of constant color can be interpreted 
as a temperature profile. 

The IHl configuration stands out from the other RB and IH models we have dis¬ 
cussed in that there is a stably stratified thermal boundary layer. The configuration thus 
provides a simple instance of penetrative convection, wherein buoyancy forces in an 
unstably stratified region drive motions that penetrate into a stably stratified region. The 
temperature field of figure reflects the dissimilarity between the unstably stratified 
upper boundary layer and the stably stratified lower one. Mixing of the cold upper layer 
with the warmer interior is accomplished by buoyantly driven cold plumes. At large R, 
the cold lower layer also can mix with the warmer interior. This mixing is driven by 
shear forces, rather than by buoyancy, and it occurs when the interior turbulence pulls 
cold eddies off the bottom boundary layer. 


3.2.2 Maximum temperatures, mean temperatures, 
and asymmetry 

The mean fluid temperature, 5(7’), behaves in IHl convection much like it does in IH3 
convection, assuming its maximum value in the static state and falling as R is raised. 
On the other hand, the mean temperature change between the boundaries, 5T, differs 
completely between the two configurations. Whereas in IH3 5T behaves rather like 
5(r), in IHl it is identically zero. The role that 5T plays in IH3 is instead approximated 
in IHl by Tmax, the maximum value that T{z) assumes over the layer. Whereas 5T 
equals the mean upward conduction across t he entire layer, Tmax captures the mean 
outward conduction, as described in ) 1.6.4.1 The quantity Tmax has been reported in 


many studies of IHl since it is easier to estimate in the laboratory than 5(7’). However, 
Tmax does not arise as easily as 5(7’) in analytical expressions. 

In addition to 5(7’) and Tmaxi IHl convection is naturally characterized by the ex¬ 
tent to which the flow creates asymmetry between upward and downward heat fluxes. 
This asymmetry can be simply conveyed by the mean fractions of produced heat leav- 
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Figure 3.2: (a) Mean temperature profiles, T (z), from the 3D DNS of Goluskin & van 
der Poel 150) for a fluid with Pr=l between no-slip boundaries. The rightmost prohle 
is that of the static state. The others, from right to left, are for R = 10^, 10^, 10*, 10^, 
and lO'**. (b) An interferogram from the experiments of Kulacki and Goldstein l84) 
with R = 1.5 • 10^ and Pr = 5.8. Any curve of constant color acts approximately as a 
graph of instantaneous temperature averaged over one horizontal direction. 
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ing across the top or bottom boundaries—or ^b- As described in ^ 1.6.2 these 
fractions are related to the dimensionless convective flux, (wT), by 


= \ + (wT) 


= (wT). 


(3.5) 


One can speak equivalently in terms of (wT), ^t, or ^b- Here we focus on ^b 
because it comes the closest to having a power-law dependence on R as in the regimes 
studied. 

Despite their simplicity and analytical attractiveness, neither 5{T) nor ^b has re¬ 
ceived much attention, although both quantities have been mentioned. Most authors 
have instead spoken in terms of top and bottom Nusselt numbers, here called Nj and 
Nb- The most common dehnitions of these numbers are 


Nj '.= 



T max 



(3.6) 


The above expressions are not normalized to be unity in the static state; instead, both 
are equal to 4. Data on Nj and Nb are typically ht to powers of R. 

To keep measures of asymmetry and temperature as separate as possible, we prefer 
not to examine Nj and Nb- Instead, we use ,^b as a measure of asymmetry and use Tmax 
and 5{T )—or their inverses, N and N —as measures of temperature. One undesirable 
feature of Nb is that it can initially drop below its static value as R is raised since 
,^B can initially fall faster than Tmax- For instance, this occurs in the data of Kulacki 
& Goldstein ll84ll . Such behavior prevents Nb from being well ht by a power of R 
near onset, and it is unlike the behavior of the RB Nusselt number, which cannot be 
smaller than its static value. Another disadvantage of using Nt and Nb is that their 
R-dependence differs only in regimes where ,^b is changing signihcantly. If the decay 
of J^B stops, as in the 2D simulations of Goluskin & Spiegel ET], then Nt and Nb will 
both be dominated by the scaling of 1 jTmax, slight changes in the asymmetry will 
not be captured well. 

We would like to summarize past data on R-dependence with hts of the form 

Tmax-^aR-^ 5{T)r^bR-l^ (3.7) 


Fits of the above form have been reported for all three quantities in ll5Qll and for 5{T) 
in EH- In two other studies where the original data are available to us ll84lll54L we 
have calculated fits to Tmax and ,^b- For the remaining studies, only hts to Nj and Nb 
are available. In these cases, we use the relations 


T 


max 


1 

Nt +Nb 


■^B = 


Nb 

Nt +Nb 


(3.8) 


The reported power-law hts to Nt and Nb dehne curves for Tmax and ,^b that are not 
pure pow ers o f R, so we have re-ht pure power laws to the latter curves. 

Table 3.2 summarizes power-law hts to the R-dependence of Tmax, 5(T), and ,^b- 
The hts to Tmax are also stated in terms of N and Ra, and the hts to 5{T) are also stated 
in terms of N and Ra- For the IHl conhguration, we have dehned these diagnostic 
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computed directly from data for a few studies I84JI154J|5T]I501 . while for other studies we have computed them from reported fits to Nt 
and Nb (see text). Simulations employ no-slip boundary conditions, except when specified otherwise. 
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quantities in 


§|TA4][r^as 


A^=l/8 7’^^ 
7V= 1/125(7’) 


Ra = R/N (3.9) 

^a = R/N. (3.10) 


The fits (3.7 1 to T^ax and 5{T) imply fits to and N and N of the form (3.4i, where 
7 = a/(l - a), 5 = i3/(l -/3), c = and d = (125)-i/(i-/^). 


3.2.2.1 Maximum temperatures 


Fits of the form T„ 


aR “ are shown in Table 


3.2 In all laboratory experiments 


and all simulations with no-slip boundaries, the exponent a lies between 0.180 and 
0.224. This means that the dimensional maximum temperature, Tmax grows with the 
volumetric heating at rates between and //0.820 ^^jg fQj. which both 

Tmax and 5(7’) are reported ll50l . the decay of Tmax is slightly faster than the decay of 
5(7’), the fit exponents being a — 0.217 and j3 = 0.204, respectively. This makes sense 
since Tmax initially must ‘catch up’ to 5(7’) as the temperature profile flattens. When 
the Tmax fits are restated in the form N ~ cRa^, the exponents range from y = 0.220 to 
7 = 0.289. The bottom end of this range is smaller than any exponents found for the 
ordinary RB1 Rayleigh number, except at very small Pr \ 


3.22.2 Mean temperatures 

The quantity 5(7’) has been reported only in two numerical studies, and each gives a fit 
of the form 5(7’) ^bR^^ for Pr = 1. Despite one study being 3D and the other 2D, the 
growth rates of 5{T) with R are very similar, having exponents of j3 = 0.204 in 3D ll50l 
and p = 0.200 in 2D ifSTIl . This is reminiscent of the Nusselt number in RB convection, 
which is not much affected by dimensionality unless Pr is small 111 15111051 . 

The dimensional mean temperature, 5{T)A, grows with the volumetric heating pro¬ 
portionally to //f* in 3D and to //O.soo jjj 2D. When the 5{T) fits are restated in the 

_ g 

form N ^ dRa , the exponents are 5 = 0.256 in 3D and 5 = 0.250 in 2D. These 5 
values are within the range of Nusselt number growth rates seen in RB 1 convection, 
though they are at the lower end of that range (cf. §3.1.2| l. We cannot yet draw com¬ 
parison with IH3 convection, for which no data on 5{T) have been reported. 

3.2.2.3 Asymmetry 

The asymmetry between upward and downward heat fluxes in IHl convection, as quan¬ 
tified by the fraction of heat that flows downward, .^b, seems to have no analogues in 
our other five IH or RB configurations. First, this fraction changes with R much more 
slowly than any other integral quantity we have discussed. Second, the P-dependence 
of .^B can differ greatly between 2D and 3D, even when Pr is not small. This is be¬ 
cause shear, rather than buoyancy, is the mechanism responsible for mixing the cooler 
lower boundary layer with the warmer interior. The asymmetry is generally greater in 
3D than in 2D because the shear-driven mixing, which helps heat escape across the 
bottom boundary, is less effective in 3D ifSOl . 
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A particularly simply question without an obvious answer is: as R ^ what is 
the limit of The extreme possibilities of either 0 or 1/2 seem most likely, although 
intermediate values are also plausible. In the highest-/? simulation data available in 
3D, falls monotonically as R is raised Il50ll . The highest-/? data available in 2D are 
quite different, except perhaps at large Pr ED. For instance, in the 2D simulations of 
Goluskin & Spiegel ED with Pr = I, the fraction reaches a minimum of 0.33 near 
R = 10^ and then increases as R is raised further. This non-monotonic /?-dependence in 
2D is yet another way that .^b stands apart from other quantities we have considered. 

When .^B decreases monotonically as R is raised, as in all past 3D studies and some 
2D ones, we can seek hts of the form .^b ~ eR^^. Table 3_2 summarizes these fits, all 
of whose decay rates are quite small. The decay exponents range from e = 0.0494 
to e = 0.0988. It remains a mystery whether such decay will continue or reverse at 
larger R 

The dependence of .^b on Pr has been examined in two studies EUEQ). The 
value of .^B seems to fall monotonically as Pr is raised, meaning that the asymmetry 
increases, until saturating at large Pr. The effect of Pr on the asymmetry is fairly 
strong—stronger than its effect on 5(7’) or Tmax- 


3.2.3 Scaling arguments 

Several scaling arguments have been put forth to explain the parameter-dependence of 
mean temperatures in IH convection ll24l l25l l26l fSD . In RB convection, the Nusselt 
number displays a wide diversity of scaling behavior in different regions of parameter 
space 11301. It is likely that the same is true of mean temperatures in IH convection 
since (inverses of) these temperatures have many parallels to the RB Nusselt number. 
This remains to be conhrmed by a wider exploration of parameter space. If such a 
diversity of scaling can indeed be found in IH convection, then any broadly applicable 
scaling arguments must reflect this. For the standard RBI configuration, the only ar¬ 
guments that attempt to capture the full range of scaling behavior are those put forth 
by Grossman & Lohse in ll52l and subsequent papers (see 11301 1. The arguments of 
15^ carry through analogously for IH convection ED- When the predicted scalings are 
phrased in terms of N and Ra, or N and Ra, they are the same as the scalings predicted 
for the Nusselt number in the RBI case. However, further work on scaling arguments 
is perhaps premature until data are available across a wider swath of parameter space. 


3.3 Future directions 


In the future study of IH convection, the main task accessible to mathematical analysis 
is proving parameter-dependent bounds on key integral quantities. The only results 
of this kind are the /?-dependent lower bounds on volume-averaged temperatures de¬ 
scribed in ) 2.3 We have conjectured in ) 1.6.3| that there should also exist /?-dependent 
upper bounds on the mean convective flux, {wT). These would amount to lower bounds 
on the fraction of heat flowing downward in IHl and on the mean temperature differ¬ 
ence between the boundaries in IH2 and IH3. Bounds are lacking also for the maximum 
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horizontally averaged temperature, Tmaxi that has often been measured in IHl exper¬ 
iments. Bounds depending analytically on the Prandtl number are highly desirable as 
well. 

There is much fertile ground for physical and computational experiments on IH 
convection. This is especially true for computation since most prior results are sev¬ 
eral decades old, so modem computers would be able to probe unexplored parameter 
regimes with relative ease. Neither the IH2 nor IH3 configuration has been simulated 
in 3D, and the DNS carried out in 2D has not approached the large R that are now 
computationally accessible. The IHl configuration has been the subject of two DNS 
studies in 3D Ill54ll50l . but a much wider exploration of parameter space is called for. 
The asymmetry between upward and downward heat fluxes in IHl is particularly hard 
to predict; even its value as R approaches infinity is not certain. In each of the three 
IH configurations, simulating a wide range of R and Pr would produce a more global 
picture of bow key integral quantities depend on the control parameters. The compli¬ 
cated parameter-dependence of Nusselt numbers in RB convection |[T] 1 13011 suggests 
that fitting integral quantities to pure powers of R will not suffice. 

A combination of mathematical analysis, simulation, and physical experimenta¬ 
tion will lead to a better understanding of the three internally heated configurations we 
have studied in this SpringerBrief. We hope that this, in turn, will lead to a better un¬ 
derstanding of more complicated occurrences of IH convection. The many past studies 
of RB convection should prove useful in guiding future studies of IH convection, and 
to this end we have described a number of analogies between the two classes of flows. 
Still, the analogies are not perfect, and some consequences of internal heating cannot 
be foreseen. Judging by the complexity of RB convection, we expect that these novel 
aspects of IH convection will remain rich areas of inquiry for many more years. 
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